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Abstract—A recent Atlantic-wide 
tag-recapture experiment run by 
the International Commission for the 
Conservation of Atlantic Tunas was an 
opportunity to directly validate otolith 
increment deposition rates for bigeye 
tuna (Thunnus obesus) and yellowfin 
tuna (T: albacares) in the region. Age 
and time at liberty were estimated 
by using annual and daily increment 
counts for sectioned otoliths from sam- 
pled fish previously injected with oxy- 
tetracycline and later recaptured. The 
use of annual increment counts resulted 
in greater age estimates than those from 
daily increment counts for fish >55 cm 
straight fork length (SFL). Use of daily 
increment counts led to underestima- 
tion of time at liberty for fish >55 cm 
SFL at recovery, compared with known 
times at liberty. In contrast, predictions 
based on annual increment counts are 
accurate across the entire size range 
of sampled fish, validating the notion 
that increments are deposited annually. 
We therefore recommend that counting 
annual increments be the preferred 
method for aging yellowfin and bigeye 
tuna from the Atlantic Ocean and that 
the use of daily increments for aging be 
limited to young of the year. Aging fish 
accurately is important for stock assess- 
ments in which data on age and growth 
play an increasingly essential role in 
examining population dynamics. It is 
crucial that otolith reading practices 
and analyses based on age data reflect 
the most up-to-date recommendations 
for age estimation. 
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Fish otoliths (or ear stones) are met- 
abolically inert structures whose con- 
centric growth rings have been widely 
used as indicators of fish age. Although 
the exact physiological and biochemical 
processes of increment formation in oto- 
liths are not well understood, it is gener- 
ally accepted that the rate of deposition 
is regulated by variations in both biotic 
(e.g., growth, feeding, reproduction, 
and stress) and abiotic (e.g., light and 
water temperature) factors (Morales- 
Nin, 2000). For many species, the otolith 
increment deposition rate has been vali- 
dated. However, for tropical bigeye tuna 
(Thunnus obesus) and yellowfin tuna 
(T. albacares), validation studies have 
been limited (Table 1), and the protocols 
accepted for age determination are not 
the same across stocks and are applied 
inconsistently across the size spectrums 
of stocks. 

The results of many studies indi- 
cate that the technique of using daily 
increment counts is not suitable for 
species that are medium- to long-lived. 
The outer daily increments become 
very narrow once fish reach a certain 


size, making them difficult to discern 
(Campana, 1992; Jones, 1992). In addi- 
tion, otolith increments may not be 
deposited daily after fish reach a cer- 
tain age or size (Neilson and Campana, 
2008). Annual increments, on the other 
hand, are generally visible through- 
out the life of a fish; therefore, counts 
of these increments are often deemed 
more accurate than daily increment 
counts, especially once a fish reaches 
a certain age or size (Casselman, 1983; 
Francis et al., 1992; Williams et al., 
2013). However, until recently, sci- 
entists thought that fish inhabiting 
tropical environments could not be 
aged reliably by using this technique 
because their environment lacks clear 
seasonal signals. Increased expertise 
and results of recent studies validating 
the rate of deposition of annual incre- 
ments in otoliths of fish in tropical 
environments have shifted that notion 
(Newman et al., 1996; Cappo et al., 
2000; Pilling et al., 2000; Morales-Nin 
and Panfili, 2005; Fowler, 2009). 
Annual growth increments in oto- 
lith sections have successfully been 


Table 1 
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Summary of studies, adapted from Williams et al. (2013), in which the periodicity of the formation of growth increments in otoliths 
was directly validated for yellowfin tuna (YFT) (Thunnus albacares) and bigeye tuna (BET) (T: obesus) sampled in the Atlantic, 
Pacific, and Indian Oceans. The age validation methods used in these studies to determine the deposition rate of daily and annual 
increments in otoliths included oxytetracycline (OTC) mark-recapture experiments, captive experiments (captive), bomb radio- 
carbon dating (““C), and strontium chloride (SrCl,) mark-recapture experiments. LM=light microscope; SEM=scanning electron 
microscope; SFL=straight fork length. 


Region 
(ocean) Species 
Eastern 
Pacific 
Eastern 
Pacific 
Western 
Pacific 
Western 
Pacific 
Eastern 


53 


74 


12 


3 


83 


method 
OTC 
OTC 

Captive 
OTC 


OTC 


Sample Validation Increment 
size 


type 

Daily 
Daily 
Daily 
Daily 


Daily 


Time at 
liberty 
Reading method (d) 


Whole otolith, LM 3-389 


Whole otolith, LM <515 


Whole otolith, LM 3-39 
Sectioned otolith,SEM 21-175 
and LM 


Sectioned otolith, LM 10-412 


Length 
range 


Age 


range 


(cm SFL) (years) Source 


40-110 


<148 


25—40 


39-91 


44—95 


Wild and Foreman, 
1980 
Wild et al., 1995 


Yamanaka, 1990# 
Lehodey and 


Leroy® 
Hallier et al., 2005 


Atlantic 

North 2 
Pacific 

Eastern 70 
Pacific 

Western 10 
Pacific 

Western 
Indian 

Western 
Indian 

Western 12 =e 
Atlantic 

Western 34 =e 
Atlantic 


Captive Daily 


OTC Daily 
SrCl, Annual 
OTC Daily 
OTC Daily 

Annual 


Annual 


Whole otolith, LM 

Sectioned otolith, LM 
Sectioned otolith, SEM 207-2420 
Sectioned otolith, LM 
Sectioned otolith, LM 
Sectioned otolith, LM 


Sectioned otolith, LM 


24-30 52 Uchiyama and 
Struhsaker, 1981 

Schaefer and 
Fuller, 2006 


Farley et al., 2006 


15-551 = =38—-135 


79-159 


3—-1166 46-142 Sardenne et al., 
2015 
Sardenne et al., 
2015 
Andrews et al., 
2020 
Andrews et al., 


2020 


8-969 48-135 


130-175 3-17 


100-180 2-18 


* Yamanaka, K. L. 1990. Age, growth and spawning of yellowfin tuna in the southern Philippines. Indo-Pac. Dev. Manag. 
Program., IPTP Work. Pap. 21, 87 p. 


used to age a number of tuna and billfish species (Gunn 
et al., 2008; Griffiths et al., 2010; Farley et al., 2013; 
Wells et al., 2013; Secor et al., 2014; Farley et al.’; Lang 
et al., 2017; Murua et al., 2017; Farley et al.?; Pacicco 
et al., 2021). The annual rate of increment deposition has 
been recently validated for bigeye and yellowfin tuna in 
the Atlantic Ocean by using bomb radiocarbon dating 
(Table 1) (Andrews et al., 2020; Pacicco et al., 2021) and 
for bigeye tuna in the western Pacific Ocean by using 
strontium chloride in mark-recapture experiments 


Farley, J., N. Clear, D. Kolody, K. Krusic-Golub, P. Eveson, and 
J. Young. 2016. Determination of swordfish growth and matu- 
rity relevant to the southwest Pacific stock. West. Cent. Pac. 
Fish. Comm. WCPFC-SC12-2016/SA-WP-11, 90 p. [Available 
from website. ] 

2 Farley, J., K. Krusic-Golub, P. Eveson, N. Clear, F Roupsard, 
C. Sanchez, 8. Nicol, and J. Hampton. 2020. Age and growth of yel- 
lowfin and bigeye tuna in the western and central Pacific Ocean 
from otoliths. West. Cent. Pac. Fish. Comm. WCPFC-SC16-2020/ 
SA-WP-02, 27 p. [Available from website.] 


(Farley et al.°). Direct comparison of ages based on 
annual increment counts versus daily increment counts 
have also been carried out for tropical tuna species, and 
the results indicate that age estimates from daily incre- 
ment counts are negatively biased compared with age 
estimates from annuli (Griffiths et al., 2010; Williams 
et al. 2043): 

In the eastern Pacific Ocean, daily increment counts are 
used for aging yellowfin and bigeye tuna up to ~150 cm 
in straight fork length (SFL) (Wild and Foreman, 1980; 
Schaefer and Fuller, 2006; Minte-Vera et al., 2020; Xu 
et al., 2020). No direct aging is carried out for fish larger 
than 150 cm SFL, even though the longevity for bigeye 
tuna has been estimated to be at least 15-16 years on the 
basis of tagging data (Langley et al., 2008). 


3 Farley, J., K. Krusic-Golub, N. Clear, P. Eveson, N. Smith, and 
J. Hampton. 2019. Project 94: workshop on yellowfin and bigeye 
age and growth. West. Cent. Pac. Fish. Comm. WCPFC-SC15-2019/ 
SA-WP-02, 14 p. [Available from website.] 
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In the western and central Pacific Ocean, Indian Ocean, 
and Atlantic Ocean, the ages of yellowfin and bigeye tuna 
are now routinely estimated by using a combination of 
counting daily increments, hereafter referred to as daily 
aging (restricted to small fish <80 cm SFL or at age 1), 
and counting annual increments, hereafter referred to as 
annual aging (Farley et al.*; Farley et al.”; Allman et al., 
2020; Pacicco et al., 2021). 

Differences in aging methods result in divergent assump- 
tions regarding the life history of fish of the same species. 
The oldest aged specimens of yellowfin tuna from the 
Atlantic Ocean and from the western and central Pacific 
Ocean have been estimated to be 18 years old (Andrews 
et al., 2020) and 15 years old (Farley et al.”), respectively, 
and both ages are based on annual aging. In contrast, the 
oldest aged specimen of yellowfin tuna from the eastern 
Pacific Ocean is estimated to be 4 years old, an age based 
on daily aging (Wild, 1986). Although differences in max- 
imum age could result from geographical differences in 
historical fishing pressure (i.e., age truncation of exploited 
populations caused by sustained, size-selective fishing), 
in this case, they are more likely the result of differences 
in aging protocols (1.e., using annual versus daily incre- 
ment counts). With fish aging underpinning estimates of 
age composition, growth, and natural mortality in stock 
assessments, it is particularly important that precise and 
unbiased aging protocols be used across all stocks. 

Age validation and comparison studies for tuna species 
are limited in the Atlantic Ocean. Our goals in this study 
were to provide additional evidence for the periodicity of 
daily and annual increment formation in yellowfin and 
bigeye tuna in the Atlantic Ocean and to provide guidance 
regarding the utility of daily and annual increment counts 
for studies of the age and growth of tropical tuna. We ana- 
lyzed a subset of otoliths sourced from a large-scale tropical 
tuna tagging campaign that was carried out from 26 June 
2015 through 28 February 2021. As part of the Atlantic 
Ocean Tropical tuna Tagging Program (AOTTP) of the Inter- 
national Commission for the Conservation of Atlantic Tunas 
(ICCAT) (AOTTP Coordination Team”), over 9000 tropical 
tuna were injected with oxytetracycline (OTC), a chemical 
marker commonly used to validate the rate of deposition of 
daily and annual increments in otoliths (Campana, 2001). 

Once a fish is injected with OTC, the chemical is quickly 
incorporated into the otolith structure and leaves a per- 
manent mark on the increment that formed at the time 
of tagging. Upon recovery, the mark can be detected by 
using fluorescence microscopy. Comparing the number of 
increments present beyond the OTC mark with the known 


* Farley, J., K. Krusic-Golub, P. Eveson, N. Clear, P. L. Luque, 
I. Artetxe-Arrate, I. Fraile, I. Zudaire, A. Vidot, R. Govinden, 
et al. 2021. Estimating the age and growth of bigeye tuna (Thun- 
nus obesus) in the Indian Ocean from counts of daily and annual 
increments in otoliths. Indian Ocean Tuna Comm. IOTC-2021- 
WPTT23-18_Rev1, 28 p. [Available from website.] 

° AOTTP Coordination Team. 2021. ICCAT Atlantic Ocean Tropical 
tuna Tagging Programme (AOTTP)—final narrative report, 57 p. 
Int. Comm. Conserv. Atl. Tunas Secr., Madrid, Spain. [Available 
from website.] 


time at liberty of a fish allows one to validate the rate of 
deposition of individual increments or growth bands. Our 
specific objectives therefore were as follows: 1) to test the 
frequency of deposition of daily and annual increments in 
the sagittal otoliths of yellowfin and bigeye tuna marked 
with OTC, 2) to compare age estimates based on daily 
versus annual increment counts for otolith thin sections 
taken from the same fish that were OTC marked and not, 
and 3) to test whether the sectioning plane used influ- 
enced total counts of daily increments and counts of incre- 
ments after the OTC mark. 


Material and methods 
Otolith sampling 


As part of the AOTTP, 3146 yellowfin tuna and 1967 big- 
eye tuna were injected with OTC. Of those injected fish, 
498 yellowfin tuna and 384 bigeye tuna have been phys- 
ically recovered to date (Suppl. Fig. 1) (online only). The 
OTC-marked fish were measured and dissected to obtain 
biological data (e.g., length, weight, and sex), and hard 
parts were extracted, cleaned, and stored for further anal- 
ysis. A subsample of fish that included the most valuable 
samples (i.e., the largest fish and fish with times at liberty 
sufficient for increments to be detected beyond the OTC 
mark) were selected for the age analyses (Fig. 1). The sub- 
sample comprised 31 bigeye tuna and 38 yellowfin tuna 
marked or recaptured across the tropical and subtropi- 
cal Atlantic Ocean, including waters off Brazil, Azores, 
Canary Islands, West Africa, St. Helena, and South Africa. 
Additional otoliths from large (127-172 cm SFL) yellowfin 
tuna not marked with OTC were donated to the project, 
2 otoliths by the Centre for Environment Fisheries and 
Aquaculture Science (St. Helena, UK) and 11 otoliths by 
the University of Cape Town in South Africa (Fig. 1). 


Otolith preparation 


Otoliths were imaged prior to sectioning and weighed to 
the nearest milligram if unbroken. The most complete oto- 
lith from each pair was selected for preparation, and the 
core was marked prior to it being embedded in Polyplex 
Clear Ortho Casting Resin (Allnex, Frankfurt, Germany). 
Otoliths were then set in resin blocks, oriented to allow 
a transverse section to be cut from the center of the oto- 
lith (Suppl. Fig. 2A) (online only). Sectioning otoliths on the 
transverse plane allowed results and methods to be directly 
compared to those of previous studies on yellowfin tuna in 
the Atlantic Ocean (Stéquert et al., 1996; Shuford et al., 
2007) and the western and central Pacific Ocean (Farley 
et al.”) and on yellowfin and bigeye tuna in the Indian 
Ocean (Stéquert and Conand, 2000; Sardenne et al., 2015). 
It also enabled us to conduct annual and daily aging on the 
same otolith because, once annual increment counts and 
OTC-mark examination were completed, the otolith could 
be ground thinner for daily increment analysis, leaving the 
remaining otolith to be retained for other purposes. 


50 100 
SFL at recovery (cm) 


250 500 750 1000 


Time at liberty (d) 
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50 100 
SFL at recovery (cm) 


250 500 750 1000 
Time at liberty (d) 


Figure 1 


Frequency distributions for straight fork lengths (SFLs) at time of recovery and for times at lib- 
erty of yellowfin tuna (YFT) (Thunnus albacares) and bigeye tuna (BET) (T. obesus) sampled in 
the tropical and subtropical Atlantic Ocean from 2015 through 2021. Most fish (31 bigeye tuna and 
38 yellowfin tuna) were injected with oxytetracycline (OTC) and recovered as part of the Atlantic 
Ocean Tropical tuna Tagging Program. An additional 11 otoliths from bigeye tuna and 2 otoliths 
from yellowfin tuna (127-172 cm SFL) were not marked with OTC. 


One section approximately 500 pm thick was cut 
through the center of each embedded otolith, ensuring 
that the primordia and core area remained within the 
otolith section. Cuts were made by using an IsoMet 1000 
precision cutting machine (Buehler Ltd., Lake Bluff, IL) 
and a single diamond wafering blade (102 x 0.3 mm) set 
to a cutting speed of 7.5 Hz (450 rpm). To allow detec- 
tion of OTC marks and annual aging, sections were 
mounted on microscope slides (76.2 x 25.4 mm) by using 
thermoplastic resin (Cystalbond 509, Electron Micros- 
copy Sciences, Hatfield, PA), with the side of the section 
farthest away from the core facing up. Each section was 
ground to a thickness between 320 and 350 pm with wet- 
and-dry sandpaper (800 and 1200 grit), lubricated with 
distilled water. Once annual increment counts and the 
OTC-mark examination were completed, the transverse 
sections were further ground down, by using 1200-grit 


wet-and-dry sandpaper, to a thickness of approximately 
60—90 ym to reveal the daily increment structure. The 
sections were polished with 2-um aluminum oxide slurry 
against a felt pad, rinsed, and dried. 

For the donated otoliths sampled from large fish, slides 
were prepared only for use in annual aging for the follow- 
ing reasons: the discrepancy between annual and daily 
increment counts already reported for otoliths from large 
specimens caught in the Pacific Ocean (Williams et al., 
2013); the value of preserving prepared slides for annual 
age estimation training purposes, especially considering 
that otoliths from very large fish are difficult to source; 
and the need to ensure that the otoliths remaining from 
each pair are still available for future analyses or valida- 
tion work (e.g., age estimation through near-infrared spec- 
troscopy or age validation through the use of bomb 
radiocarbon dating). 
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Age estimation with otoliths 


Annual increment counts The methods used for the annual 
aging of bigeye tuna and yellowfin tuna follow those 
developed for other tuna species that are routinely aged 
at Fish Ageing Services, in Queenscliff, Australia (Farley 
et al., 2006; Gunn et al., 2008; Farley et al., 2013), and 
those in information available in the literature at the time 
(i.e., Lang et al., 2017). Information on aging bluefin tuna 
sampled in the Atlantic Ocean was also used as a basis 
to aid interpretation of what may constitute an annual 
growth zone in bigeye tuna and yellowfin tuna from the 
Atlantic Ocean (Rodriguez-Marin et al., 2007; Neilson and 
Campana, 2008; Rodriguez-Marin et al., 2014; Secor et al., 
2014). Before reading each section, the ground surface 
of the otolith was covered in a thin layer of low-viscosity 
immersion oil (Type A, Cargille-Sacher Laboratories Inc., 
Cedar Grove, NJ) to fill in any residual scratches and aid 
in the imaging process. Sections were examined at 25x 
magnification under transmitted light with a Leica M125 
C routine laboratory stereo microscope (Leica Microsys- 
tems, Wetzlar, Germany). 

The annual age of each section was estimated by count- 
ing opaque growth zones (which appear dark under trans- 
mitted light). The last fully completed opaque zone before 
the otolith edge was counted only if translucent otolith 
material was detected between the outer edge of the last 
opaque zone and the otolith edge. All age readings were 
made without knowledge of fish size, otolith weight, sex, 
location of capture, or time at liberty. A single TIFF image 
was captured, and the distances between the primordium 
and the outer edge of each opaque zone was measured. For 
each section, a readability score of 1-5 was assigned, with 
1 meaning a zone pattern could not be interpreted and 5 
meaning there was a clear pattern of alternating opaque 
and translucent zones. Also recorded were any relevant 
comments related to the otolith structure or interpreta- 
tion of the zones. 


Daily increment counts Daily increments were counted 
by using a Leitz Diaplan compound microscope (Leica 
Microsystems) with transmitted light at various magni- 
fications ranging between 400x and 1000x, depending 
on the area of the otolith being interpreted. The inter- 
pretation of the otolith microstructure in tuna can be 
subjective, and it can often be difficult to distinguish 
sub-daily increments from the assumed daily incre- 
ments. These difficulties were recognized by Shuford 
et al. (2007) and Sardenne et al. (2015), and, as per their 
methods, the sub-daily increments were characterized 
by faint and incomplete rings. Further complicating the 
daily increment count was the presence of subsections 
within some otoliths for which little or no increment pat- 
tern could be detected. Because daily increment widths 
tend to be autocorrelated (Campana, 1992), the counts of 
micro-increments within these subsections were interpo- 
lated from the surrounding areas. The method used for 
the interpretation of daily increments is consistent with 
those methods published for reading transverse sections 


(Lehodey and Leroy®; Shuford et al., 2007; Sardenne 
et al., 2015) and frontal sections (Schaefer and Fuller, 
2006). The otolith sections were read at least 2 times by 
the single reader before a final count was completed and 
the number of daily increments was recorded. 


Validation of increment deposition rate: 
oxytetracycline-mark detection 


Each otolith section from a fish that was injected with 
OTC was examined for the presence and position of the 
OTC mark by using a Leitz Diaplan compound microscope 
fitted with a 100-W incident ultraviolet light source and a 
Leitz D filter block (Leica Microsystems, with excitation 
filter of 450-520 nm) to suit the fluorescent properties of 
the OTC. Images were taken at magnifications of 25x and 
40x, and, if the OTC mark was very faint, at 100x magni- 
fication with a DFK 31AU03 digital camera (The Imag- 
ing Source, Charlotte, NC) attached to the microscope and 
the camera’s corresponding image analysis software (IC 
Measure, vers. 2.0.0.245; The Imaging Source). Scale bars 
appropriate to the magnification were included on each of 
the images. Images were collected for 2 reasons: 1) to cap- 
ture an image of the section when the OTC mark was at 
its brightest, given that continual exposure to UV and the 
further processing of sections for daily increment counts 
can result in the OTC mark becoming very faint or even 
disappearing within the section, and 2) to detect, through 
the use of measured distances, the approximate position of 
the OTC mark in the otolith section, when not viewing the 
section under ultraviolet light. To make this detection, the 
distances between the OTC mark and the otolith margin 
on both the inside and outside of the ventral arm were 
measured by using the image capture software. 

The number of annual increments observed between 
the position of the OTC mark (based on the measured 
distances taken in the previous step) and the otolith 
edge were recorded following the annual aging protocols 
detailed in the previous section. With annual aging, the 
direct comparison between the count of annual increments 
after the OTC mark and the time at liberty can be mis- 
leading because, unlike with daily increment counts, age 
estimates based on counts of opaque zones do not provide 
a fractional age (e.g., a fish at liberty for only 4 months 
could have an opaque zone between the OTC mark and 
the edge if the opaque zone completed formation shortly 
after tagging). It is possible to estimate a fractional age if 
the timing of annual zone formation is well known. How- 
ever, for yellowfin and bigeye tuna sampled in the Atlantic 
Ocean, annual aging is only in its infancy, and the timing of 
zone formation for fish caught in the equatorial and south- 
east Atlantic Ocean is currently unknown. Therefore, to 


® Lehodey, P., and B. Leroy. 1999. Age and growth of yellowfin 
tuna (Thunnus albacares) from the western and central Pacific 
Ocean as indicated by daily growth increments and tagging 
data. Standing Comm. Tuna Billfish 12, Work. Pap. YFT-2, 
21 p. Ocean. Fish. Progr., Secr. Pac. Community, Noumea, New 
Caledonia. [Available from website.] 


provide an alternative, yet practical, comparison between 
estimated time at liberty based on annual aging and the 
true time at liberty, we calculated an adjusted estimate 
(in fractional years) for each of the OTC-marked sections 
using methods modified from those used by Cappo et al. 
(2000). This adjusted count was based on the following: 
1) the position of the first opaque zone after the OTC 
mark, relative to the OTC mark, 2) the number of annual 
zones observed after the mark, and 3) the relative dis- 
tance between the last assumed annual opaque zone and 
the otolith edge—assuming that one translucent and one 
opaque zone are formed each year. 

The sections prepared for daily increment counts were 
reexamined, and increment counts were made from the 
OTC mark to the otolith edge along a count path similar 
to that used for the total daily increment counts. For the 
sections in which the OTC mark had become undetectable 
because of additional otolith processing, the positions of 
the OTC marks were determined from the measurements 
taken on each of the otolith sections during the imaging 
procedure. 


Longitudinal sections 


To examine whether the sectioning plane that was used 
influenced increment counts, longitudinal (frontal) sections 
of the remaining otolith for a subset of fish (6 bigeye tuna 
and 5 yellowfin tuna) were prepared following methods 
used by Schaefer and Fuller (2006) (Suppl. Fig. 2B) (online 
only). These longitudinal sections were used to age the sub- 
set of fish following methods used by Williams et al. (2013), 
in which opaque zones were counted from the primordia 
to the otolith margin along the clearest count path in the 
area of the section adjacent to the proximal edge. The total 
count of daily increments and the number of daily incre- 
ments between the OTC mark and the edge were recorded. 


Analytical methods 


For the comparison of age estimates, daily and annual 
increment counts from the same individuals were plotted 
against each other to help visualize differences. For the 
exercise of validating the increment deposition rate, daily 
increment, annual increment, and adjusted annual incre- 
ment counts after the OTC mark were plotted against 
time at liberty to aid visualization of differences between 
estimated and true times at liberty. Age-agreement tables 
and age-bias plots (Ogle, 2016) were constructed, and an 
Evans—Hoenig test of symmetry (Evans and Hoenig, 1998) 
was applied to evaluate bias in the times at liberty esti- 
mated by using daily increment and adjusted annual incre- 
ment counts. The Evans—Hoenig test is designed to detect 
bias in paired-age data and, in our study, was applied to 
paired data for estimated and true times at liberty. A sig- 
nificant P-value (<0.05) indicates that differences between 
estimated times at liberty and true times at liberty were 
due to bias and not random error, implying that daily and 
annual increments were not strictly deposited on a daily 
and annual basis. For these calculations, times at liberty 


Fishery Bulletin 121(1—2) 


were converted to integer month to reduce the number of 
categories being compared. 

To determine whether the sectioning plane that was 
used influenced total and post-OTC-mark counts of daily 
increments, the increment counts from transversely sec- 
tioned otoliths were compared with counts from fron- 
tal sections of otoliths. A Bland—Altman plot (Bland and 
Altman, 1999) was used to evaluate bias and define the 
interval of agreement between reads (total counts) from 
otoliths sections made with the 2 cutting planes. 


Results 
Age estimation: daily versus annual increment counts 


Detailed results pertaining to each otolith analyzed 
for daily and annual increment counts are presented in 
Supplementary Tables 1 and 2 (online only). The annual 
increments in the form of one opaque and one translucent 
zone, although often difficult to interpret from sub-annual 
marks (also known as check marks or split zones), were 
apparent in all transverse sections examined, except in 
sections of otoliths from certain young-of-the-year samples 
for which the first zone had not yet formed. Zone struc- 
ture was observed on both the dorsal and ventral side of 
otolith sections, although our experience in aging otoliths 
from samples of other tuna species indicates that the zone 
structure within the ventral arm is generally easier to 
interpret in comparison to that within the dorsal arm and 
that counts of assumed annual increments on the ventral 
side are likely to provide for more accurate age estimates. 
Annual increment counts ranged from 0 to 8 for bigeye tuna 
and from 0 to 10 for yellowfin tuna. A transverse section of 
an otolith from one of the oldest yellowfin tuna examined 
in this study (not OTC marked) appears in Figure 2, which 
shows the positions of each of the opaque zones used to 
assign age based on counts of annual increments. 

Clear daily increments, consisting of alternating opaque 
and translucent zones, were observed throughout most of 
the ventral arms of the otolith sections examined. How- 
ever, for some areas in many of the otolith sections exam- 
ined, interpretation of the daily increment pattern was 
subjective. Overlapping and split increments often made 
interpretation difficult, and no increments were visible in 
some regions within the otolith sections. The number of 
regions that were difficult to interpret within an otolith 
generally increased with an increase in relative otolith 
size. Similar to that observed in otoliths from yellowfin 
tuna caught in the Pacific Ocean (Lehodey and Leroy®), the 
core area in otoliths examined in this study consisted of a 
primordium and approximately 8—10 fine daily increments 
(with widths of ~2 ym) followed by a check mark. Incre- 
ment widths became progressively wider, and at approxi- 
mately the fifteenth increment, growth increments became 
very optically dense and their thickness increased up to 
20 um at their maximum concavity. Counts from the first 
obvious zone to the first apex averaged between 30-35 
increments for both species. The internal zones closer to 
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Figure 2 

Image of a section of an otolith from the oldest individual aged in 
this study: an Atlantic Ocean yellowfin tuna (Thunnus albacares) 
captured on 18 October 2018, measuring 161 cm straight fork length 
and weighing 80.3 kg. The white circles indicate the positions of the 
opaque zones used to assign age based on counting of annual incre- 
ments. The letter N indicates the position of the nucleus. The opaque 
zone on the otolith margin (black arrow) is incomplete and was not 
included in estimation of integer age. This otolith was donated to 
the study by the Centre for Environment Fisheries and Aquaculture 
Science and was not marked with oxytetracycline. 


the sulcus were more regularly spaced and did not have as 
much overlapping, splitting, or merging of zones compared 
with the zones close to the distal and ventral edges. Given 
these observations, it is our recommendation that the ini- 
tial count path should include the clearer internal incre- 
ments rather than run close to the distal and ventral edges 
of the otolith. 

Unfortunately, for bigeye and yellowfin tuna, after 
approximately 150-180 daily increments, the internal 
structure of otoliths becomes difficult to interpret and 
a count path closer to the outer edge of the ventral arm 
needs to be used. As per the method of Lehodey and Leroy’°, 
the reading of the outer increments was directed along the 
axis of maximum concavity of increments. In many of the 
otoliths examined, subsections along the preferred read- 
ing path had a zone structure that was either difficult to 
interpret or not present. These areas were relatively com- 
mon in otoliths from individuals of both species and often 
occurred with a change in the otolith growth plane, and 
in visual comparisons made with the otolith images cap- 
tured during preparation for aging, the positions of these 
areas often corresponded to the positions of the opaque 
zones marked on images. In these cases, the zone struc- 
ture adjacent to these areas usually had a clearer pattern 
of alternating opaque and translucent zones. For these 
cases, daily increments were counted in the adjacent 


area until increments could again be interpreted 
clearly along the preferred aging path. If the 
adjacent area was also unclear, interpolation was 
used (Suppl. Fig. 3) (online only). Daily increments 
were detectable close to the outer edge even in 
the largest otoliths from individuals of both spe- 
cies (114 and 159 cm SFL for bigeye and yellow- 
fin tuna, respectively). These increments were at 
least 1.5—2 ym wide and well above the minimum 
resolution of light microscopy. 

Age estimates, both daily increment counts 
and raw annual increment counts, are presented 
in Figure 3 and in Supplementary Tables 1 and 2 
(online only). In the count comparison exercise, 
sizes of yellowfin tuna at recapture ranged from 
45 to 159 cm SFL, daily increment counts ranged 
from 247 to 1168, and annual increment counts 
ranged from 0 to 8. For bigeye tuna, recapture 
sizes ranged from 50 to 114 cm SFL, daily incre- 
ment counts ranged from 248 to 642, and annual 
increment counts ranged from 0 to 3. Age esti- 
mates of bigeye tuna with no annual zones were 
strongly negatively biased compared with the 
corresponding daily increment counts (Fig. 3). 
However, the edge type for all 6 otolith sections 
with no annual increments was classified as 
opaque, indicating that the fish was likely closer 
to age 1 than to age 0. Although it is difficult to 
objectively compare raw annual increment counts 
to fractional ages obtained from daily increment 
counts, for fish at sizes beyond 55 cm SFL (~age 1), 
age estimates based on annual aging tended to 
be higher than estimates based on daily aging for 
both species, with differences generally more pronounced 
in fish greater than 100 cm SFL for which age estimates 
from annual aging were consistently higher than esti- 
mates from daily aging (Fig. 3). The most drastic differ- 
ence between age estimates was observed for a yellowfin 
tuna of 159 cm SFL: it had a daily increment count of 753 
and an annual increment count of 8. 


Validation with oxytetracycline marking 


Oxytetracycline marks were detected in 32 of the 38 oto- 
liths from yellowfin tuna and 20 of the 31 otoliths from 
bigeye tuna that were injected as part of the AOTTP and 
included in this study. For most examined otoliths, the 
OTC mark was clear and easy to detect from the back- 
ground fluorescence of the otolith section (Fig. 4). Lengths 
at recapture of the fish successfully marked with OTC 
ranged from 45 to 159 cm SFL for yellowfin tuna and from 
50 to 114 cm SFL for bigeye tuna (although the length 
ranges for fish with OTC marks on their otoliths are the 
same as those for all recaptured fish, not all recaptured 
fish had an OTC mark on their otolith). Known times at 
liberty ranged from 7 to 995 d for yellowfin tuna and from 
18 to 879 d for bigeye tuna. 

The estimated times at liberty (.e., counts of increments 
after the OTC mark) were plotted against the known 
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Figure 3 


Comparison of age estimates based on daily (squares) and 
raw annual (triangles) increment counts on the same oto- 
lith, plotted by straight fork length (SFL), for yellowfin 
tuna (YFT) (Thunnus albacares) and bigeye tuna (BET) 
(T. obesus) sampled in the tropical and subtropical Atlan- 
tic Ocean between 2015 and 2021. Circles in the top panel 
indicate the increment counts for the 11 YFT for which 
only annual increment counts were taken. 


times at liberty (Figs. 5 and 6). If the increments were 
truly deposited on a daily basis, and assuming there was 
no interpretation error, the points in each plot should fall 
directly on the solid black line that represents the 1:1 ratio 
of increment counts to times at liberty (Figs. 5A and 6A). 
Except for the readings of otoliths from one bigeye tuna 
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and one yellowfin tuna, the daily increment counts fell 
below the 1:1 line, indicating that counts were lower than 
the true times at liberty. Results also indicate that the 
larger the fish size, the greater the discrepancy between 
the age assumed based on daily increment count and the 
true time at liberty. The raw counts of annual increments 
indicate a high level of correspondence between the num- 
ber of opaque zones observed after the OTC mark and the 
actual time at liberty (Figs. 5B and 6B). The otolith from 
one yellowfin tuna at liberty for just over 2 years had only 
one visible opaque zone after the OTC mark, but the oto- 
lith had a wide translucent edge, indicating that a second 
opaque zone was due for deposition on the margin (Fig. 6). 
When the reader estimated times at liberty on the basis 
of the position of the OTC mark, the number of annual 
zones observed after the mark, and the relative distance 
between the last assumed annual opaque zone and the oto- 
lith margin, the estimated time at liberty agreed with the 
1:1 line more closely and no directional bias was apparent 
(Figs. 5C and 6C). 

Age-bias plots (Figs. 7 and 8) and results of the Evans— 
Hoenig test confirm our observation that daily increment 
counts are biased in that they are low in comparison to 
true times at liberty (Suppl. Tables 3-6) (online only). With 
the test, we detected a significant difference between the 
daily increment counts and the known time at liberty in 
both yellowfin tuna (y?=32, df=11, P=0.0008) and bigeye 
tuna (y7=18, df=10, P=0.05), indicating that increments 
were not systematically deposited on a daily basis for 
either species. The results from the Evans—Hoenig test 
of symmetry between the times at liberty estimated from 
annual increment counts and the actual times at liberty 
indicate that there was no systematic bias in times at 
liberty estimated from annual increment counts for both 
yellowfin tuna (y’=3.33, df=6, P=0.77) and bigeye tuna 
(y7=5.81, df=5, P=0.33). 


Transverse versus longitudinal sections 


When directly compared, the counts of increments in 
transverse versus longitudinal sections of otoliths did not 
differ greatly, nor did they indicate an obvious bias due to 
the preparation method. The Bland—Altman plot indicates 
a high degree of agreement between the counts made with 
the 2 reading methods, and there was no systematic differ- 
ence between the paired counts (Fig. 9). There was a small 
negative bias in mean difference in counts with the 2 
methods (-13 increments), a difference that is negligible 
in the context of age estimation (Fig. 9). All data points 
that were plotted fell within the 95% limits of agreement 
(-83, +57). This range indicates how far apart paired 
counts were most likely to be for most individuals, and the 
range is considered acceptable given the inter- and intra- 
reader variability typically observed in aging of tropical 
tuna. We did note, however, that daily increments within 
the longitudinal sections were easier to interpret than 
those within the transverse sections and that this differ- 
ence in readability was particularly true for sections of 
otoliths from bigeye tuna. 
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Figure 4 


Examples of otolith sections from (A and B) a bigeye tuna (Thunnus obesus) and (C and D) a 
yellowfin tuna (T: albacares) sampled in the Atlantic Ocean. The images in panels A and C were 
illuminated with transmitted light, and the yellow line indicates the approximate position of the 
oxytetracycline (OTC) mark on each section. In panels B and D, the OTC marks are clearly vis- 
ible in the corresponding sections illuminated with ultraviolet light. The white circles indicate 
the positions of the annual opaque zones used for aging. The bigeye tuna was recaptured on 29 
April 2000 after 879 d at liberty, and its lengths at release and recapture were 45 and 114 cm 
straight fork length (SFL), respectively. The times at liberty estimated for this fish from daily and 
annual increment counts were 494 and 821 d, respectively. The yellowfin tuna was recaptured on 
2 February 2020 after 995 d at liberty, and its lengths at release and recapture were 77 and 154 cm 
SFL, respectively. The times at liberty estimated for this fish from daily and annual increment 


counts were 357 and 912 d, respectively. 


Discussion 


In this study, reading otolith thin sections with transmit- 
ted light proved to be an appropriate method for obtain- 
ing age estimates based on annual increment counts for 
bigeye and yellowfin tuna sampled in the Atlantic Ocean. 
The age-at-length data presented in Figure 3 are similar 
to those data reported from other studies of tropical tunas 
that also are based on counts of annual increments over 
the age range examined, namely studies on yellowfin tuna 
in the Atlantic Ocean (Lang et al., 2017; Pacicco et al., 
2021) and on both bigeye tuna and yellowfin tuna in the 
western Pacific Ocean (Farley et al., 2006; Farley et al.”). 
The strong relationship between the estimated time at lib- 
erty and the true time at liberty and the lack of systematic 
bias confirm that counting annual growth zones on oto- 
liths from individuals of these 2 species is a valid method 
for age estimation. Furthermore, there were no instances 
of multiple annual rings being detected in a single year, 
as have been detected for little tunny (Euthynnus alletter- 
atus), skipjack tuna (Katsuwonus pelamis), and blackfin 


tuna (T: atlanticus) from the Atlantic Ocean (Adams and 
Kerstetter, 2014). 

Further work is needed to convert zone counts to deci- 
mal age estimates for both bigeye tuna and yellowfin tuna 
from the Atlantic Ocean and to increase the accuracy of 
age estimates for our geographically disparate fish sam- 
ples. Multiple methods are currently used for assigning 
fractional ages for fish species. A common approach is to 
adjust the raw count on the basis of knowledge of the tim- 
ing of band deposition and the birth date of the fish. How- 
ever, neither the timing of band deposition nor the natal 
origin of the fish sampled were known for the otoliths used 
in this analysis; therefore, this approach could not be used 
reliably. An alternative and arguably preferable approach 
developed by Farley et al.” does not require information on 
spawning and timing of band formation. Instead, this 
approach involves developing a relationship between the 
age based on daily increment counts and the length of the 
ventral otolith arm for age 0+ fish and determining esti- 
mates of average annual increment width for each age 
class. Still, this approach requires large numbers of fish to 
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Figure 5 


Comparison of (A) daily, (B) annual, and (C) adjusted annual increment counts, plotted with true 
times at liberty, for bigeye tuna (Thunnus obesus) marked and recaptured in the tropical and 
subtropical Atlantic Ocean during 2015-2021. Points that appear below the solid black line that 
represents the 1:1 ratio of increment counts to times at liberty indicate that the increment count is 
an underestimation of the true time at liberty, and points that appear above the line indicate that 
the increment count is an overestimation of the true time at liberty. Adjusted annual increment 
counts were determined on the basis of the position of the oxytetracycline mark, the number of 
annual zones observed after the mark, and the relative distance between the last assumed annual 
opaque zone and the otolith edge. Symbols relate to the quality of the recovery date: date is known 
only approximately and could be off by a few days, date is known exactly, and the quality of the 
date is unknown. Shades of gray used for symbols correspond to straight fork lengths (SFLs) of 
the fish at recapture according to the scale provided, with darker shades indicating longer lengths. 


be sampled and to be analyzed a priori, something that 
could not be accomplished within the scope of our study 
but that could certainly be attempted in the future by 
using the remaining sampled fish from the AOTTP. 

In the otoliths prepared for daily aging, areas with read- 
ily identifiable daily increments were observed throughout 
the otolith section in both transverse and longitudinal sec- 
tions, albeit they were separated by areas in which the 
increment structure was difficult to interpret or was not 
present at all. Although the theory that the otoliths of 


older fish contain daily increments that are below the res- 
olution limits of light microscopy is commonly held, in our 
study, daily increments were still detected close to the 
outer edge of otoliths even in the sections of the largest 
otoliths of both species examined. Schaefer and Fuller 
(2006) made the same observation for otoliths from bigeye 
tuna that were up to 145 cm SFL and were sampled from 
the eastern Pacific Ocean, although they used longitudinal 
sections. Similarly, Williams et al. (2013), comparing both 
longitudinal and transverse sections, detected daily 
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Figure 6 


Comparison of (A) daily, (B) annual, and (C) adjusted annual increment counts, plotted with true 
times at liberty, for yellowfin tuna (Thunnus albacares) marked and recaptured in the tropical and 
subtropical Atlantic Ocean during 2015-2021. Points that appear below the solid black line that 
represents the 1:1 ratio of increment counts to times at liberty indicate that the increment count is 
an underestimation of the true time at liberty, and points that appear above the line indicate that 
the increment count is an overestimation of the true time at liberty. Adjusted annual increment 
counts were determined on the basis of the position of the oxytetracycline mark, the number of 
annual zones observed after the mark, and the relative distance between the last assumed annual 
opaque zone and the otolith edge. Symbols relate to the quality of the recovery date: date is known 
only approximately and could be off by a few days, date is known exactly, and the quality of the 
date is unknown. Shades of gray used for symbols correspond to straight fork lengths (SFLs) of 
the fish at recapture according to the scale provided, with darker shades indicating longer lengths. 


increments in the outer part of the otolith section of the 
largest (175 cm SFL) southern bluefin tuna (7! maccoyii) 
examined in their study (senior author, personal observ.), 
despite the considerable divergence of the age estimates 
based on total annual (29 years) and daily (4 years) incre- 
ment counts. 

The advantage of using longitudinal sections is that this 
sectioning plane provides the longest axis for counting 
daily increments, with increment spacing that is usually 
wider than that in corresponding transverse sections. This 


wider spacing can lead to fewer areas within the otolith 
section having daily increments that are hard to interpret 
or requiring interpolation. Although a detailed comparison 
of sectioning methods was not an objective of this study, we 
did compare the results between a small number of oto- 
liths prepared by using the 2 different sectioning planes. 
This comparison was important to avoid potential bias in 
the results from using either of the otolith preparation 
methods and to provide insight into which sectioning 
method should be preferred for use in future aging work 
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Figure 7 


Age-bias plots showing the differences between annual and daily increment counts, or estimated times at liberty, 
against the true times at liberty for bigeye tuna (Thunnus obesus) marked and recaptured in the tropical and 
subtropical Atlantic Ocean from 2015 through 2021. In each panel, a histogram of the differences between the 
estimated and true times at liberty are provided in the right margin. The number of samples (7) for each time 
at liberty (in months) is provided in a histogram in the top margin. Circles indicate mean differences between 
the estimated and true times at liberty at each true month at liberty, and the bars stretch from the minimum to 
the maximum difference observed. The open circle indicates that the mean is significantly different from zero. 


based on counts of daily increments for these 2 species. 
Preliminary results indicate no clear bias in the daily 
increment counts between the preparation methods used 
and are consistent with those of an earlier study on yellow- 
fin tuna in the Indian Ocean for which preparation method 
also did not influence the age estimates (Stéquert et al., 
1996). We note, however, that daily increments within the 
longitudinal sections in our study were easier to interpret 
than those within the transverse sections, and this differ- 
ence in ease of detecting increments was particularly true 
for bigeye tuna. We therefore recommend that a more thor- 
ough comparison of the 2 methods be completed. 

Results from both the comparison of annual versus daily 
increment counts and the validation of deposition rates of 
daily increments in otoliths based on examination of OTC- 
marked fish indicate that counting daily increments in 
otoliths of yellowfin and bigeye tuna may lead to under- 
estimating age for fish larger than 55 cm SFL (greater 
than ~age 1). The difference between age estimates based 
on the 2 methods of counting can become quite large as 
fish grow older, as was observed for one of the yellowfin 
tuna examined in our study, with a difference in age of 6 


years. Williams et al. (2013) counted both daily and annual 
increments on their largest southern bluefin tuna and 
found an overall difference of 20 years between the age 
estimates from the 2 counts. Sardenne et al. (2015) simi- 
larly found that daily increment counts in otoliths of large 
(>100 cm FL) bigeye tuna from the Indian Ocean resulted 
in underestimates of age. 

One explanation for this underestimation could be that 
daily increments are not systematically deposited at a 
daily rate. It has been hypothesized that increment depo- 
sition may cease to be daily after a certain age or life his- 
tory stage and that changes in environmental conditions 
might disrupt the regular pattern of daily increment depo- 
sition (e.g., Francis et al. 1992). Another potential explana- 
tion relates to the fact that there are areas within otoliths 
(even in those from relatively small fish, <75 cm SFL) 
where daily increments are either extremely difficult to 
interpret or do not appear to form at all; such areas also 
have been found in other studies (Shuford et al., 2007; 
Sardenne et al., 2015). Daily increments tend to have a 
merging or splitting structure in some areas, making it 
difficult to distinguish individual increments, and other 
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Figure 8 


Age-bias plots showing the differences between the annual and daily increment counts, or estimated times at lib- 
erty, against the true times at liberty for yellowfin tuna (Thunnus albacares) marked and recaptured in the tropi- 
cal and subtropical Atlantic Ocean from 2015 through 2021. In each panel, a histogram of the differences between 
the estimated and true times at liberty are provided in the right margin. The number of samples (7) for each time 
at liberty (in months) is provided in a histogram in the top margin. Circles indicate mean differences between 
the estimated and true times at liberty at each true month at liberty, and the bars stretch from the minimum 
to the maximum difference observed. Open circles indicate that the means are significantly different from zero. 


areas are devoid of increments, forcing analysts to resort 
to interpolation to fill the missing patterns. Campana 
(1992) suggested that the use of interpolation is consid- 
ered reasonable in cases in which the number of interpo- 
lated daily increments is relatively small in comparison to 
the total increment count. Unfortunately, for the otoliths 
of bigeye and yellowfin tuna examined in this study, these 
areas were relatively common and likely the main contrib- 
uting factor in age underestimation. 

Despite these limitations to the use of daily increment 
counts for aging yellowfin and bigeye tuna, daily incre- 
ments do hold considerable value in the overall aging 
process. Daily increment counts are still important for 
confirming the location of the first annual increment in 
sectioned otoliths, as has been reported for yellowfin tuna 
in the Atlantic Ocean (Lang et al., 2017) and both bigeye 
and yellowfin tuna in the Pacific Ocean (Farley et al., 
2006; Farley et al.”), and for providing length-at-age data 
during the first year of life. More recently, aging through 
the use of daily increments counts, particularly out to a 
count of 365 d, has proven to be a valuable tool in refining 
counts of annual zones and, therefore, in providing more 


accurate fractional ages and has also allowed improve- 
ment in the growth curves of yellowfin and bigeye tuna 
from the western and central Pacific Ocean (Farley et al.”). 
With this in mind, the challenge is determining the point 
at which age estimation based on daily increment counts 
provides a larger source of error than the use of annual 
zone counts. Finding that point should be considered a 
vital step in developing suitable aging methods for any 
species, let alone for tropical tunas, and the need to obtain 
that knowledge highlights the value of conducting a large- 
scale tag-recapture program such as the AOTTP. 


Conclusions 


On the basis of our findings, the use of annual increment 
counts is the best method for aging yellowfin and bigeye 
tuna sampled in the Atlantic Ocean, and the use of daily 
increment counts should be limited to aging young of the 
year. Additional work is needed to resolve the timing of 
opaque zone formation in otoliths from fish from differ- 
ent parts of the Atlantic Ocean and to derive fractional 
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Figure 9 


Bland—Altman plot of the paired readings of total daily 
increment counts for otoliths from bigeye tuna (Thunnus 
obesus) (circles) and yellowfin tuna (T: albacares) (trian- 
gles) marked and recaptured in the tropical and subtropi- 
cal Atlantic Ocean between 2015 and 2021. The horizontal 
dashed line indicates the mean bias in the difference in 
counts between the transverse and longitudinal sections. 
The horizontal dotted lines indicate the 95% limits of 
agreement, calculated as the mean plus and minus 1.96 
standard deviations (SDs). 


ages, and we believe the AOTTP otolith collection could 
serve that purpose. Although the AOTTP has now ended, 
networks across the Atlantic Ocean are being maintained 
to ensure that there is a continued effort to recover OTC- 
marked fish and to analyze hard parts. As this effort con- 
tinues, more valuable and informative samples (i.e., larger 
fish, with longer times at liberty) will become available, 
and analysis of their otoliths is likely to confirm our find- 
ings and recommendations for these 2 key tuna species. 

Validating age estimation methods is particularly valu- 
able for improving stock assessments; the importance of 
obtaining correct age and growth information in assessing 
stock status has been well documented (Maunder and 
Punt, 2013). Lang et al. (2017) and Pacicco et al. (2021) 
have suggested that yellowfin and bigeye tuna from the 
Atlantic Ocean can be aged by using counts of annual 
increments and that the maximum longevity for both spe- 
cies was far greater than originally determined. More 
importantly, the aging protocols used in those recent stud- 
ies have been verified by the use of bomb radiocarbon anal- 
ysis (Andrews et al., 2020). The results of that recent work, 
along with the results from our analysis of otoliths from 
OTC-marked fish in our study, support the use of annual 
increments for aging yellowfin and bigeye tuna sampled in 
the Atlantic Ocean. 
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The findings of our study are particularly important 
because the recent shift to base age estimates on annual 
increment counts has led to raising the estimated maxi- 
mum age in the latest assessment of the Atlantic stock of 
yellowfin tuna from 11 to 18 years (ICCAT, 2020), which 
in turn has caused substantial changes in stock status 
(ICCAT, 2019). The ICCAT has assessed the Atlantic stocks 
of both the bigeye tuna and yellowfin tuna with modeling 
approaches of varying complexities, from simple produc- 
tion models to integrated statistical assessment models, 
and the capacity to accommodate age data in several dif- 
ferent ways. Although age data have not yet been used to 
their full extent for tuna stocks in the Atlantic Ocean, we 
anticipate that that will change in the future, given the 
recent advances in the study of age and growth of tropical 
tuna species. 


Resumen 


Un experimento reciente de marcado-recaptura en todo 
el Atlantico llevado a cabo por la Comisié6n Internacional 
para la Conservacion del Attin del Atlantico fue una opor- 
tunidad para validar directamente las tasas de incremen- 
tos de otolitos para el patudo (Thunnus obesus) y el atun 
de aleta amarilla (T. albacares) en la regi6n. Se estimaron 
la edad y el tiempo en libertad utilizando conteos de incre- 
mentos anuales y diarios de otolitos seccionados de peces 
muestreados previamente inyectados con oxitetraciclina y 
recapturados posteriormente. Los conteos de incrementos 
anuales dieron lugar a mayores estimaciones de edad que 
las realizadas con conteos diarios, para los peces mayores 
de 55 cm de longitud furcal recta (SFL). El conteo de incre- 
mentos diarios condujo a una subestimacidén del tiempo 
en libertad para los peces mayores de 55 cm de SFL al ser 
recapturados, en comparacion con los tiempos en libertad 
conocidos. Por el contrario, las predicciones basadas en los 
recuentos de incrementos anuales son precisas en todo el 
rango de tallas de los peces muestreados, validando asi, 
que los incrementos se depositan anualmente. Por lo tanto, 
recomendamos que el conteo de incrementos anuales sea 
el método preferido para determinar la edad del attin de 
aleta amarilla y el patudo del Océano Atlantico y que el 
uso de incrementos diarios para determinar la edad se 
limite a los juveniles del ano. La determinacién precisa 
de la edad en peces es importante para la evaluacion de 
poblaciones, en las que los datos sobre edad y crecimiento 
desempenan un papel cada vez mas esencial en el estudio 
de dindamica de poblaciones. Es crucial que las practicas 
de lectura de otolitos y los analisis basados en los datos de 
edad reflejen las recomendaciones mas actualizadas para 
la estimacion de la edad. 
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Abstract—Assessments of the status 
and ecological linkages of target fish 
species related to marine habitats are 
playing important roles in enhancing 
the effectiveness of marine protected 
areas (MPAs) and of management in 
other coastal areas. In this study, we 
assessed juvenile yield and the abun- 
dance of goldlined spinefoot (Siganus 
guttatus) by using local knowledge 
and visual census in coastal waters of 
Vietnam. Genetic structure was exam- 
ined by using DNA analysis of juve- 
niles, collected in seagrass beds and 
mangroves at major lagoons and estu- 
aries, and of adults, collected from coral 
reefs. More than 38.6 million juveniles 
were collected from 8 coastal lagoons 
and estuaries. The coastal lagoons and 
estuaries with seagrass beds supported 
high yield of juveniles, and the nearby 
MPAs had higher densities of adults. 
Three distinct genetic groups were 
identified, with one group found off the 
central part of the eastern coast from 
Con Co Island to Quy Nhon Bay, a sec- 
ond group found off the southeastern 
coast from Nha Trang Bay to the Con 
Dao Islands, which are located off the 
southern tip of Vietnam, and a third 
group found off the southwestern coast 
at the Phu Quoc Islands in the Gulf of 
Thailand. 
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Results from previous research indicate 
that more than 210 species of 34 reef 
fish families have pelagic stages, that 
larvae disperse away from their spawn- 
ing area (Green et al., 2015), and that 
many of them migrate from their set- 
tlement habitats of mangroves and 
seagrass beds in coastal lagoons and 
estuaries (nursery grounds) to adja- 
cent coral reefs for feeding or shelter 
(Luo et al., 2009; Honda et al., 2013; 
Le et al., 2019; Vrdoljak et al., 2021). 
In an increasing number of studies of 
marine fish species, strong genetic dif- 
ferentiation among populations from 
different regions, separated by dis- 
tances from tens to a few hundred kilo- 
meters, has been observed (Hellberg 
et al., 2002; Jones et al., 2005; Laikre 
et al., 2005; Planes et al., 2009). These 
patterns may be influenced by several 
ecological and behavioral factors, such 
as the timing of reproduction (Selkoe 
et al., 2006), larval dispersal ability 
(Kinlan and Gaines, 2003), homing 
behavior (Almany et al., 2007), ocean- 
ographic patterns (White et al., 2010), 


biogeographic barriers (McMillan 
and Palumbi, 1995), local adaptation 
(Conover et al., 2006), habitat special- 
ization (Knudsen et al., 2006), and the 
ecological requirements of each species 
(Rocha et al., 2002). 

Reports from several studies indi- 
cate that many reef fish families 
(Lutjanidae, Haemulidae, Lethrinidae, 
Scaridae, Siganidae, and others) uti- 
lize seagrass beds and mangroves in 
nearshore and estuarine regions as 
their nursery grounds (Verweij et al., 
2008). Therefore, understanding the 
degree of connectivity between popu- 
lations in different geographic areas, 
or the specific location of genetic 
breaks, sets the scale at which man- 
agement strategies for species need 
to be applied. This information is 
important in the management of reef 
fish resources in many coastal areas, 
especially in marine protected areas 
(MPAs) (Cowen and Sponaugle, 2009). 

Results of studies on the genetic 
structure of several species of siganids 
indicate that there are differences in 


connectivity among species. In waters of Japan, Iwamoto 
et al. (2009) reported a higher degree of interpopulation 
gene flow for bluntnosed spinefish (Siganus spinus) than 
for goldlined spinefoot (S. guttatus) between Okinawa 
Island, Miyako Island, and Ishigaki Island. Magsino and 
Juinio-Menez (2008) found higher genetic structuring 
in dusky spinefoot (S. fuscescens) than in streamlined 
spinefoot (S. argenteus) along the eastern coasts of the 
Philippines, and they also observed genetic structure for 
the these 2 species in the South China Sea and South 
Philippine Sea but no structure for them in the Sulu Sea 
and Inland Sea. In the South China Sea, high connectiv- 
ity between populations of reef fish species, including the 
false pennant bullfish (Heniochus acuminatus), the six- 
barred wrasse (Thalassoma hardwicke), and the threespot 
dascyllus (Dascyllus trimaculatus), has been reported 
from genetic studies, indicating population differences in 
some areas (Ablan et al., 2002; Chen et al., 2004). 

In the coastal waters of Vietnam, results of a genetic 
study indicate broad connectivity of the goldstripe sardi- 
nella (Sardinella gibbosa), the giant tiger prawn (Penaeus 
monodon), the threespot dascyllus, the cupped oyster Cras- 
sostrea rivularis, and the giant clams Tridacna crocea and 
T. squamosa between a northern region (from the Cat Ba 
Islands to Da Nang Bay), a central section (from Quy Nhon 
Bay to Nha Trang Bay), an area in and around the Con 
Dao Islands, located off the southern tip of Vietnam, and 
a region in and around the Phu Quoc Islands, located off 
the southwestern coast in the Gulf of Thailand, although 
a weak separation of populations has been found between 
the Con Dao Islands and Phu Quoc Islands (Dang et al., 
2014). Similar patterns have been reported for the ornate 
spiny lobster (Panulirus ornatus) in Southeast Asia (Dao 
et al., 2015) and for a complex of mullet species, Mugil 
cephalus, in Vietnam (Tran et al., 2017). 

The goldlined spinefoot, also known as the orange-spotted 
spinefoot and the golden rabbitfish, is 1 of 28 species in the 
Indo-Pacific family Siganidae. It is widely distributed 
in shallow-water habitats Gncluding lagoons and river 
mouths, mangrove swamps, seagrass beds, coral reefs, 
and sandy and rocky bottoms) to a depth of 6 m in the 
eastern Indian Ocean and the Indo-West Pacific Ocean 
(Woodland, 1990; Duray, 1998; Froese and Pauly, 2021). 
This species is one of the most economically important 
commercial fishes in tropical and subtropical regions. It 
has attracted attention from the aquaculture industry 
in the Indo-Pacific and in the Indian Ocean (Lam, 1974; 
Randall et al., 1990). Juveniles settle in seagrass beds in 
lagoons and estuaries during certain periods of the year 
(Kami and Ikehara, 1976), whereas adults are found in 
coral reefs (Froese and Pauly, 2021). Because this species 
has a complex life cycle with linkages to different habitats 
at different stages, understanding the genetic structure 
of populations of the goldlined spinefoot is important for 
developing effective strategies for its proper and sustain- 
able management. 

On the western side of the South China Sea, the coast- 
line of Vietnam spans 14° of latitude and extends for 
3260 km. The coastal waters support more than 156,000 ha 
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of mangroves (Do et al., 2005) and 17,000 ha of seagrass 
beds (Luong et al., 2012), and coral reefs are mainly dis- 
tributed along the coast and around the adjacent islands 
over a total area of about 13,400 ha (Nguyen and Vo, 2014; 
Nguyen et al., 2019). In Vietnam, juvenile goldlined spine- 
foot have been recruited to seagrass beds and mangrove 
swamps in the coastal lagoons and estuaries in the central 
part of the country from Thua Thien Hue Province (which 
includes Tam Giang—Cau Hai and Lap An Lagoons) to 
Binh Dinh Province (where Thi Nai Lagoon is located) 
(Le and Le, 2006). Large numbers have been observed set- 
tling during the dark phase of the moon in the summer 
months from May to August (from day 20 in one month to 
day 5 in the next month) and sometimes in to September 
(Le and Le, 2006; Nguyen and Mai, 2018). Adults have 
been found in the major coral reefs along the coast and 
around nearshore islands (Nguyen and Mai, 2020). Large 
schools of juvenile goldlined spinefoot have been caught 
by purse seine, lift net, stow net, and long trap cage 
along the coast of Vietnam, including in Thi Nai Lagoon 
(8.58 million juveniles in 2008) (Nguyen et al., 2010) and 
in Thu Bon estuary (780 kg, equivalent to 7.02 million 
juveniles, in 2015) (Nguyen and Mai, 2018). However, data 
at a larger spatial scale on the distribution of juveniles 
and broodstocks is lacking, and such information is needed 
for effective management of this resource in Vietnam. 

The aims of this study were as follows: 1) to assess the 
juvenile yield and adult abundance (broodstock) of gold- 
lined spinefoot; 2) to determine the genetic diversity, 
structure, and linkages of juvenile goldlined spinefoot col- 
lected in coastal estuaries and lagoons and of broodstock 
collected in coral reefs in coastal waters of Vietnam; and 
3) to provide a baseline for future management decisions 
regarding this species in Vietnam. 


Materials and methods 
Sampling strategies 


Assessments of juvenile yield Assessments of the yield of 
juvenile goldlined spinefoot were conducted through con- 
sultations with local communities. The word juvenile is 
used herein for the specimens of goldlined spinefoot, with 
total lengths (TLs) of 18-45 mm, collected by fishermen 
for the aquaculture industry. Although there are some 
limits to accuracy, information and data collected through 
local consultations, including local ecological knowledge, 
have been widely used in many countries around the 
world where accurate scientific data are insufficient 
(Begossi, 2015; Berkstrém et al., 2019). Data from local 
consultations have been considered semi-quantitative 
data for rapid assessments of the status of habitats and 
related resources that contribute to providing a compre- 
hensive baseline for appropriate planning and sustain- 
able utilization of resources in Latin American countries 
(Fischer et al., 2015). The local consultations were con- 
ducted at 8 fishing communities near major estuaries and 
lagoons along the coast of Vietnam in November and 
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December 2019 under a national project. The 8 locations For each consultation, biological scientists from the 
where consultations were conducted are Cua Tung estu- Institute of Oceanography, Vietnam Academy of Science 
ary, Cua Viet estuary, Tam Giang—Cau Hai Lagoon, Lap and Technology (VAST), together with local authorities 
An Lagoon, Thu Bon estuary, Tra Bong estuary, Sa Ky selected 10-20 consultants, including experienced fish- 
estuary, and Nha Phu Lagoon (Fig. 1). | ermen that use different fishing gears, local dealers, and 
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Figure 1 


A map of the locations where local consultations regarding goldlined spinefoot (Siganus guttatus) 
(black circles) or coral reef surveys of reef fish taxa (open circles) were conducted along the coast 
of Vietnam. Yield of juveniles (18—45 mm total length [TL]) was assessed through local consulta- 
tions at 8 fishing communities near major estuaries or lagoons in November and December 2019. 
Density of adult (>100 mm TL) goldlined spinefoot was assessed at 11 locations at coral reefs in 
established or proposed conservation areas during July-September 2020 as part of this study (Con 
Co Island, Hai Van—Son Cha, the Ly Son Islands, and Nha Trang Bay) or from 2017 through 2021 
under support from other projects (the Cu Lao Cham Islands, Quy Nhon Bay, Nui Chua National 
Park, Hon Cau Island, the Phu Quy Islands, the Con Dao Islands, and the Phu Quoc Islands). For 
DNA analysis, goldlined spinefoot were sampled from 2017 through 2020: juveniles from 3 loca- 
tions in an estuary or a lagoon where consultations occurred, as well as from Thi Nai Lagoon, and 
adults from 9 locations at coral reefs where surveys were conducted. The numerals in parentheses 
indicate the number of specimens that were collected at each of these 13 sites. 
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marine aquaculturists, to provide the needed information. 
Questions were focused on fishing activities for the juve- 
nile yield of goldlined spinefoot and designed to gather the 
following information: fishing gear and seasons; number 
of fishing boats; number of fishermen per fishing boat; and 
catch per unit of effort (CPUE), defined as catch per boat 
per day or night. 


Assessments of adult abundance The density of adult gold- 
lined spinefoot was assessed together with those of other 
groups of reef fishes at 11 widely distributed locations at 
coral reefs that have been established as MPAs (Con Co 
Island, the Cu Lao Cham Islands, Ly Son the Islands, Nha 
Trang Bay, Nui Chua National Park, Hon Cau Island, the 
Con Dao Islands, and the Phu Quoc Islands) or are being 
proposed as MPAs (Hai Van—Son Cha and the Phu Quy 
Islands) or locally managed marine area (Quy Nhon Bay) 
in Vietnam. The recorded sizes of adult goldlined spinefoot 
collected on coral reefs were >100 mm TL. Census surveys 
were conducted at 4 locations—3 sites off the central part 
of the eastern coast, Con Co Island, Hai Van—Son Cha, and 
the Ly Son Islands, and 1 site farther south in Nha Trang 
Bay—during July-September 2020 within the framework 
of this project and at 7 other locations—2 sites off the 
central part of the eastern coast in 2017 at the Cu Lao 
Cham Islands and in Quy Nhon Bay, 3 sites farther south 
at Nui Chua National Park in 2018, Hon Cau Island in 
2020, and the Phu Quy Islands in 2021, 1 site at the Con 
Dao Islands in 2018, and 1 site at the Phu Quoc Islands 
in 2019—under support from other projects (Table 1, 
Fig. 1). At each study site, 2 transects, each 100 m long, 
were placed haphazardly in 2 reef zones (1 transect on 
a reef flat and 1 transect on a reef slope) parallel to the 
shoreline. The depths below low tide level at transects 
were approximately 2—4 m at reef flats and 5-12 m at reef 
slopes. Each 100-m transect was divided into 4 replicates 
(each 20 m long by 5 m wide). Hence, 8 replicate transects 
were surveyed at each study site. 

Along the transect line, the fish observer used scuba gear 
and swam slowly to count the number of individuals and to 
estimate the size (TL) of each individual fish from visual 
estimates using the methods proposed by English et al. 
(1997) and Hodgson et al. (2006). Fish were identified to 
species level by using Randall et al. (1990). The survey time 
at each site was about 50-60 min. The surveys were all car- 
ried out during daylight hours between 0900 and 1400 local 
time. To avoid disturbance from divers, fish surveys were 
conducted 15 min after the first transect was placed. 


Assessments of genetic diversity and structure The sequence 
data used in this study are from fish specimens that 
were collected from 4 locations at estuaries and lagoons 
and from 9 locations at coral reefs in the coastal waters 
of Vietnam from 2017 through 2020 under support from 
different projects (Fig. 1). Three locations in an estuary 
or a lagoon, including Tam Giang—Cau Hai Lagoon, Tra 
Bong estuary, and Thi Nai Lagoon, and 4 locations at coral 
reefs, including Con Co Island, Hai Van—Son Cha, the Ly 
Son Islands, and Quy Nhon Bay, were sampled from 2019 
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through 2020 as part of this study. Thu Bon estuary and 
the Cu Lao Cham Islands were sampled in 2017 for a 
study supported by the VAST, and 4 other reef locations, 
including Nha Trang Bay, Nui Chua National Park, the 
Con Dao Islands, and the Phu Quoc Islands, were sam- 
pled from 2018 through 2019 as part of a study funded by 
the U.S. Agency for International Development. Because 
information on occurrence and distribution of juvenile 
goldlined spinefoot at Nha Phu Lagoon and other lagoons 
and estuaries off the southern coast was not yet available 
prior to the designing of sampling for this study, samples 
were not collected at these locations. 

Fish samples were collected by using cast nets and 
by purchasing specimens at fish landing sites. At each 
location, 29-33 goldlined spinefoot were collected in 
estuaries or lagoons (29-105 mm TL) and at coral reefs 
(202-340 mm TL). A total of 401 individuals (125 juveniles 
and 276 adults) were sampled. 

For each individual, a 1-cm? section of the right pectoral 
fin and a 1-cm” section of the upper caudal fin were cut out 
with fine scissors, washed with water, and immediately 
fixed in a solution of 95% alcohol in a 2-mL screw-cap tube 
labeled with species name, code, and preservation solu- 
tion. Then, these fin samples were cooled and stored fro- 
zen at —20°C for further analysis in the laboratory. 

Samples of DNA were extracted by using the cetyl- 
trimethylammonium bromide method (Adamkewicz and 
Harasewych, 1996) and preserved in Tris-EDTA buffer 
solution (10 mM Tris-HCl, pH 7.2, 1 mM EDTA, pH 8.0) 
at —20°C. The concentration and purity of the extracted 
DNA were tested on 1.5% agarose gel and on a spectropho- 
tometer (BioSpec-nano’, Shimadzu Corp., Kyoto, Japan). 
The use of microsatellites or single sequence repeats was 
referenced from existing studies of molecular biology and 
population genetics for goldlined spinefoot. In this study, 
5 microsatellites used for population genetic analyses of 
goldlined spinefoot and suggested by Mao et al. (2016) 
were applied: HLZY-3 (F: TAACGGTTCTATCAGGG and 
R: TGCTTCGGATTCAGG at 55°C), HLZY-5 (F: TTCAT- 
CACTGCCTGTCCTT and R: AGCGTGTCATTGTTG- 
GGT at 62°C), HLZY-12 (F: C@GACCTGCTTGAGT and 
R: AACCTGTTTGGTGGAA at 47°C), HLZY-13 (F: TGG- 
CAGTGATGAGGTG and R: CATTTGGAATCGCAGT at 
57°C), and HLZY-21 (F: CTGACTCCCAACTTCC and R: 
CAGACCTGTTTCCAATC at 52°C). 

Five forward and reverse primer pairs were run through 
the polymerase chain reaction (PCR) process, in which 
primers were marked by HEX or FAM fluorescence stain- 
ing at the 5’ end. The PCR amplification was performed 
with a BOECO Thermal Cycler TC-PRO (Boeckel and Co. 
GmbH and Co., Hamburg, Germany). The PCR reaction 
was performed at a total volume of 50 pL, including 10 pL 
of buffer, 0.2 pL Bioline MyTaq HS Polymerase (Meridian 
Bioscience, Cincinnati, OH), 0.4 1M of each primer, and 5 ng 
of the template DNA. The PCR process for microsatellites 


1 Mention of trade names or commercial companies is for identi- 
fication purposes only and does not imply endorsement by the 
National Marine Fisheries Service, NOAA. 
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Table 1 


The number of sites at which visual census surveys for adult goldlined spinefoot (Sig- 
anus guttatus), among other reef fish species, were conducted at each of 11 locations at 
coral reefs along the coast of Vietnam during 2017-2021. At the time of this study, all 
locations either had been established or were being proposed as marine protected areas 
(MPAs) or locally managed marine area. An asterisk (*) indicates that the location was 


an established MPA. 


No. of sites 
Location surveyed 
Con Co* 10 
Hai Van—Son Cha 10 
Cu Lao Cham* 10 
Ly Son* 10 
Quy Nhon Bay 24 
Nha Trang Bay* 10 
Nui Chua* 10 
Hon Cau* 12 
Phu Quy tt 


Time of surveys 


August 2020 
July 2020 
July 2017 

September 2020 
April 2017 
July 2020 

September 2018 

August 2020 
April 2021 


Sources 


This study 

This study 
Nguyen et al.’ 
This study 
Nguyen? 

This study 

Mai and Nguyen® 
Hoang et al.* 
Nguyen et al.° 


Con Dao* 9 
Phu Quoc* 10 


June 2018 
April 2019 


Mai and Nguyen? 
Nguyen® 


‘ Nguyen, V. L., T. H. Dao, X. D. Mai, T. K. H. Phan, and H. T. Nguyen. 2020. Studies on 
population connectivity of target fishes among marine habitats in the World Biosphere 
Reserve of Cu Lao Cham—Hoi An. Final report for the Vietnam Academy of Science 
and Technology project VAST06.02/17-18, 24 p. [In Vietnamese.] [Available from Inst. 
Oceanogr., Vietnam Acad. Sci. Tech., 01 Cau Da St., Nha Trang 650000, Vietnam.| 


2 Nguyen, V. L. 2017. Additional assessments and identification of locations with 

high marine biodiversity in the Quy Nhon Bay. Final report of the project CRSD/ 
BD/3.d.4.1/2016 under the provincial project “Resources in coastal waters for sustainable 
development (CRSD) in Binh Dinh Province,” 142 p. [In Vietnamese.] [Available from Inst. 
Oceanogr., Vietnam Acad. Sci. Tech., 01 Cau Da St., Nha Trang 650000, Vietnam.] 


3 Mai, X. D., and L. V. Nguyen. 2021. Status and trends of change of reef fish communi- 
ties at some locations in the coastal waters of south Vietnam in the period of 2017-2019. 
Technical report of the project titled, “Study on coral reef resilience in comparative areas 
in south Vietnam for marine biodiversity conservation in a changing world” (project code: 
AID-OAA-A-11-0012), 60 p. [In Vietnamese.] [Available from Dep. Mar. Living Resour., 
Inst. Oceanogr., Vietnam Acad. Sci. Tech., 01 Cau Da St., Nha Trang 650000, Vietnam.] 


* Hoang, X. B., T. K. H. Phan, K. H. Phan, X. D. Mai, and M. Q. Thai. 2020. Status 

and trends of change in marine biodiversity and related resources at key sites in the 
Hon Cau Marine Protected Area. Technical report of the project titled, “Advice on the 
planning of review, re-zoning, and boundary of Hon Cau Marine Protected Area, Binh 
Thuan Province,” 57 p. [In Vietnamese.] [Available from Dep. Mar. Living Resour., Inst. 
Oceanogr., Vietnam Acad. Sci. Tech., 01 Cau Da St., Nha Trang 650000, Vietnam. | 


° Nguyen, T. H., M. Q. Thai, X. D. Mai, and B. X. Hoang. 2022. Assessments of status of 
marine biodiversity and related resources in key marine habitats (coral reefs, seagrass 
beds) in the waters surrounding Phu Quy district. Technical report of the project titled 
“Advice on the development of the project for the establishment of Phu Quy Marine Pro- 
tected Area, Binh Thuan Province” (project code: 01-DABTBPQ) 50 p. [In Vietnamese. ] 
[Available from Dep. Mar. Living Resour., Inst. Oceanogr., Vietnam Acad. Sci. Tech., 01 
Cau Da St., Nha Trang 650000, Vietnam.] 


° Nguyen, V. L. 2019. Assessments of marine biodiversity and re-zoning for manage- 
ment of Phu Quoc Marine Protected Area. Final report of provincial project submitted 
to the Department of Agriculture and Rural Development of Kien Giang Province, 
310 p. [In Vietnamese.] [Available from Inst. Oceanogr., Vietnam Acad. Sci. Tech., 01 
Cau Da St., Nha Trang 650000, Vietnam.] 


or single sequence repeats had six stages. The first stage 
was the initial denaturation at 95°C for 1 min, followed by 
35 cycles of denaturation at 95°C for 15 s, an annealing 
phase at a suitable temperature for each primer for 15 s, 
and an extension phase at 72°C for 90 s. Then, the final 
extension phase began at 72°C for 10 s and ended at 20°C 


for 60 s. The PCR products were then run on 1.6% aga- 
rose gel/100 mL TBE buffer (0.5x) with 1.0 pL gelled and 
detected by ultraviolet transilluminator (UVP DigiDoc-It, 
Analytik Jena GmbH, Jena, Germany) prior to genotyping 
on an Applied Biosystems 3730XL DNA Analyzer (Thermo 
Fisher Scientific, Waltham, MA). 
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Data storage and analysis 


Annual yield of juveniles was calculated as follows: num- 
ber of boats multiplied by fishing time multiplied by aver- 
age yield per time unit (CPUE). 

The density of fish at each study site was calculated 
as the mean value of the 8 replicates from the 2 tran- 
sects on the reef flat and reef slope. To improve the homo- 
geneity of variances, raw data on abundance of adults 
was log, .(x+1) transformed prior to the analyses. When 
the analysis of variance detected significant differences, 
Tukey’s honestly significant difference test was applied 
in order to separate differences between locations. 

The software GeneMarker, vers. 2.7.4 (SoftGenetics, State 
College, PA), was used to determine the length of alleles. 
Null alleles were detected by using the software Micro- 
Checker, vers. 2.2.3 (Van Oosterhout et al., 2004). Genetic 
variability estimates for each locus, such as those made by 
using the Hardy—Weinberg equilibrium, number of alleles 
on each satellite, number of rare alleles, observed heterozy- 
gosity, and expected heterozygosity were calculated by using 
the software GenAlIEx 6.1 (Peakall and Smouse, 2006). The 
statistical significance tests for these indices were revised 
on the basis of the Bonferroni correction (Rice, 1989). 

Evaluation of differences between groups in genetic 
structure (fixation index [F's7] and pairwise F'g7) was done 
through analysis of molecular variance with Arlequin, 
vers. 3.11 (Excoffier et al., 2005). We used the model-based 
clustering method in the software Structure, vers. 2.3.4 
(Pritchard et al., 2000), to elucidate the group structure 
and visualized results in images by using the program 
Structure Harvester, vers.0.6.94 (Earl and vonHoldt, 
2012). In the model of the software Structure, a Bayesian 
approach is applied to estimate the probability of correctly 
assigned genotypes in the data in the number of clusters 
(K). Each value of K (1-10 clusters) was run for 10 repe- 
titions, and the average log probability of the data, L(K), 
was calculated. The L(K) plateaued (or 
continued to increase slightly) when 
K approached the true value, and high 
variance between runs also occurred 
(Rosenberg et al., 2001). The number of 
K groups was then evaluated by using 
AK, an ad hoc statistic based on the sec- 
ond order rate of change of L(K) (Evanno 
et al., 2005). 
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Juvenile yield 


Information for the local consultations 
at 8 locations in estuaries and lagoons 
indicates that juvenile goldlined spine- 
foot were collected only from Tam 
Giang—Cau Hai Lagoon (in Thua Thien 
Hue Province) to south-central Vietnam, 
with 38.61 million juveniles (equivalent 
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to 4290 kg) collected in 2019. Among the 8 locations, the 3 
most important locations, including Tam Giang Cau—Hai 
Lagoon, Thu Bon estuary, and Nha Phu Lagoon, supported 
high CPUE (215,000—826,000 juveniles per boat per day) 
and yield (9.77—13.01 million juveniles). The locations of 
Lap An Lagoon and Sa Ky estuary supported the lowest 
CPUE (0.03—0.04 juveniles per boat per day) and yield 
(<1 million juveniles at each location) (Fig. 2). In gen- 
eral, juvenile goldlined spinefoot were found in lagoons 
and estuaries that had seagrass beds, with an exception 
of Tra Bong estuary. However, there was no correlation 
between the area of seagrass beds and both CPUE and 
yield of juveniles of this species (Fig. 2). Information from 
local fishermen indicates that juvenile goldlined spinefoot 
were also found in the Cua Tung and Cua Viet estuaries 
(in Quang Tri Province). However, because they occurred 
in very low numbers, the local fishermen did not make an 
effort to collect them. 


Abundance of adult stock 


Results of an analysis of data from 138 sites at 11 reef 
locations in the coastal waters of Vietnam that were 
surveyed during 2017-2021 indicate that the overall 
mean density of goldlined spinefoot was 0.37 fish/100 m7 
(standard error of the mean [SE] 0.04), with the sites in 
MPAs supporting a significantly higher mean density 
(0.51 fish/100 m? [SE 0.06]) than that of the sites that 
were not in MPAs (0.09 fish/100 m? [SE 0.03]) (analysis 
of variance: P<0.001). Among the sites in MPAs, the mean 
densities in the Cu Lao Cham, Nha Trang Bay, and Phu 
Quoc MPAs were 3.3-51.2 times higher than and signifi- 
cantly different (analysis of variance: P<0.001) from those 
in other MPAs (the Con Co, Ly Son, and Nui Chua MPAs) 
and those in sites not in MPAs (Hai Van—Son Cha and Quy 
Nhon Bay) (Fig. 3). No goldlined spinefoot were recorded at 
the Hon Cau MPA, the Phu Quy Islands (proposed MPA), 
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Figure 2 


The relationship between catch per unit of effort (CPUE) and yield of juvenile 
goldlined spinefoot (Siganus guttatus) in major estuaries and lagoons along 
the coast of Vietnam, based on local consultations conducted in 8 fishing com- 
munities during November and December 2019. The black diamonds indicate 
the locations that had seagrass beds and the extent of these habitats. 
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Figure 3 


Mean density of adult goldlined spinefoot (Siganus guttatus) sampled at 11 locations at major 
coral reefs in the coastal waters of Vietnam between 2017 and 2021. An asterisk (*) indicates that 
a location was an established marine protected area at the time of this study. Error bars indicate 
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standard errors of the means. 


and the Con Dao MPA. The reef locations that supported 
higher densities of adult goldlined spinefoot (the Cu Lao 
Cham and Nha Trang Bay MPAs) than other reef locations 
are closer geographically to the lagoons and estuaries 
with higher juvenile yield (Thu Bon estuary and Nha Phu 
Lagoon). However, sites not in MPAs (Hai Van—Son Cha 
and Quy Nhon Bay) supported much lower densities than 
MPA sites despite their proximity to Tam Giang—Cau Hai 
and Thi Nai Lagoons, which had very high juvenile yield 
for collections made in this study (Fig. 3). 


Genetic diversity and structure 


Results of analysis of sequence data for 347 samples from 
401 specimens collected from 13 locations indicate that 
the total number of alleles of 5 microsatellites was 110, of 
which 382 alleles were private (Table 2). Among the 5 micro- 
satellites, HLZY-12 had the highest alleles per locus (36), 
followed by 25, 22, 18, and 9 alleles per locus for HLZY- 
18, HLZY-5, HLZY-21, and HLZY-3, respectively. The pop- 
ulation of goldlined spinefoot in the Phu Quoc Islands 
appears to be the most diverse, given that samples from 
specimens caught at that location had the highest num- 
ber of alleles (73) compared with the numbers of alleles 
for samples from specimens collected at the other loca- 
tions (41-57). The samples from fish caught at Phu Quoc 
Islands also had the highest number of private alleles (20). 
The values of observed heterozygosity and expected het- 
erozygosity of each single sequence repeat varied slightly 
between locations, ranging from 0.33 to 1.00 and from 0.49 
to 0.92, respectively (Fig. 4). In general, there was no devi- 
ation from the Hardy—Weinberg equilibrium in sequence 
data for most locations, with some exceptions for micro- 
satellite HLZY-3 (Thu Bon estuary, Thi Nai Lagoon, and 
the Phu Quoc Islands), HLZY-5 (Hai Van—Son Cha, the Cu 
Lao Cham Islands, Thi Nai Lagoon, Quy Nhon Bay, and 
Nha Trang Bay), HLZY-12 (the Cu Lao Cham Islands), 


HLZY-13 (Con Co Island, the Con Dao Islands, and the 
Phu Quoc Islands), and HLZY-21 (Hai Van—Son Cha, Tra 
Bong estuary, and Nui Chua National Park) (Table 3). 


Genetic linkages among marine habitats 


The results of the analyses of molecular variance indicate 
a significant difference in genetic structure (F'g7=0.015, 
P<0.05) among genetic groups, a central group (specimens 
from Con Co Island to Quy Nhon Bay), a southern group 
(specimens from Nha Trang Bay to the Con Dao Islands), 
and a southwestern group (specimens from the Phu Quoc 
Islands), at the sampling locations with a low percentage of 
genetic variation (1.5%; Table 4). Results from comparison 
of pairwise Foy values indicate that genetic differentiation 
occurred at all locations (at levels of 1.2-6.5%) and that 
the difference was significant (P<0.05). Generally, results 
from analyses of molecular variance reveal that there 
were 2 distinct groups in the genetic structure of collected 
goldlined spinefoot among the 13 sampling locations: a 
central group that included specimens from Con Co Island 
to Quy Nhon Bay and a southern group that included fish 
caught from Nha Trang Bay south and all the way around 
to the Phu Quoc Islands on the southwestern coast, with 
the genetic similarity among locations within each group 
reaching 99-100%. This genetic similarity indicates that 
genetic linkages occurred only within each group. 

The results from Bayesian analysis (done with the soft- 
ware Structure) indicate that AK ranged from 0.018305 to 
0.774668 when K was set from 1 to 8. The clearest peak in 
AK occurred when K was set to 3, with a value of 0.774668, 
indicating that 3 distinct genetic clusters (assigned the col- 
ors of green, red and blue in visualizations) existed in gold- 
lined spinefoot among the sampling locations (Fig. 5). The 
central group included specimens from Con Co Island to 
Quy Nhon Bay, with most fish having the majority of their 
genome assigned to the green cluster (Fig. 5). The southern 
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Table 2 


Total number of alelles and private alleles for each microsatellite found in sequence 
data from goldlined spinefoot (Siganus guttatus) sampled from 2017 through 2020 in 
the coastal waters of Vietnam. Sites where fish were collected include the following 
locations at estuaries, lagoons, or coral reefs: Con Co Island (CCO), Tam Giang—Cau 
Hai Lagoon (TAG), Hai Van—Son Cha (HVS), Thu Bon estuary (THB), the Cu Lao Cham 
Islands (CLC), Tra Bong estuary (TBO), the Ly Son Islands (LSO), Thi Nai Lagoon 
(TNA), Quy Nhon Bay (QNH), Nha Trang Bay (NTR), Nui Chua National Park (NCH), 
the Con Dao Islands (CDA), and the Phu Quoc Islands (PQU). 


Microsatellite 


Sampling eae er BPs SS Ce ol al oF 2 PE a en 
location HLZY-3 HEZY-5 HLZY-12 HLZY-13 HLZY-21 Total 


Alleles 
CCO 
TAG 
HVS 
THB 
CLC 
TBO 
LSO 
TNA 
QNH 
NTR 
NCH 
CDA 
PQU 
Total 

Private alleles 
THB 
CLC 
LSO 
QNH 
PQU 
Total 
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group was composed of individuals caught at loca- 
tions from Nha Trang Bay to the Con Dao Islands, 
forming the red cluster (Fig. 5). The southwestern 
group, consisting of specimens caught off the Phu 
Quoc Islands, formed the blue cluster (Fig. 5). 


Discussion 
The results of this study indicate that the central | HLZY-3 HLZY-5 HLZY-12 WHLZY-13  HLZY-21 
coast of Vietnam supported the highest yield of Microsatellite 


juvenile goldlined spinefoot in 2019, with higher 


iel ded at T Gi - Hai L 
yield secprned A aes ‘ne an, Han Lagoon, Mean values of observed heterozygosity (Ho) and expected heterozygos- 
Thu Bon estuary, Thi Nai Lagoon, and Nha Phu . 
: ity (He) for each microsatellite used in genetic analysis of goldlined spine- 
Lagoon compared to that at the other estuaries foot (Siganus guttatus) sampled from 2017 through 2020 in the coastal 
and lagoons. This greater level of yield may be waters of Vietnam. Error bars indicate standard errors of the means. 
related to the presence and structure of seagrass 


beds (Nagelkerken, 2009; Nguyen et al., 2015). In 


Figure 4 


this study, CPUE had no correlation with areas Higher densities of adult goldlined spinefoot were 
of seagrass beds; therefore, differences in the composition mainly found during coral reef surveys at the locations 
of seagrass species may contribute to marked differences that have been established as MPAs rather than at loca- 
in the total yield of juvenile goldlined spinefoot among the tions where MPAs do not exist. This finding is consistent 


lagoons and estuaries. with those of other studies in which higher densities or 
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Table 3 


P-values from tests for deviation from the Hardy—Weinberg equilibrium 
(HWE), by location for each microsatellite found in sequence data from 
goldlined spinefoot (Siganus guttatus) sampled from 2017 through 2020 
in the coastal waters of Vietnam. An asterisk (*) indicates a significant 
difference (P<0.05) from HWE. Sites where fish were collected include 
the following locations at estuaries, lagoons, or coral reefs: Con Co Island 
(CCO), Tam Giang—Cau Hai Lagoon (TAG), Hai Van—Son Cha (HVS), 
Thu Bon estuary (THB), the Cu Lao Cham Islands (CLC), Tra Bong estu- 
ary (TBO), the Ly Son Islands (LSO), Thi Nai Lagoon (TNA), Quy Nhon 
Bay (QNH), Nha Trang Bay (NTR), Nui Chua National Park (NCH), the 
Con Dao Islands (CDA), and the Phu Quoc Islands (PQU). 


Microsatellite 
Sampling 


location 


CCO 
TAG 
HVS 
THB 
CLC 
TBO 
LSO 
TNA 
QNH 
NTR 
NCH 
CDA 
PQU 


HLZY-3 


0.45 
0.59 
0.91 
0.02* 
0.27 
0.95 
0.95 
0.01* 
0.99 
O29 
0.62 
0.99 
0.00* 


HLZY-5 


0.81 
0.63 
0.03* 
1.00 
0.00* 
0.48 
1.00 
0.01* 
0.00* 
0.00* 
022 
0.73 
0.84 


HLZY-12 


0.82 
0.89 
0.60 
0.74 
0.00* 
0.49 
0.68 
0.53 
0.30 
0.45 
0.35 
0.54 
0.74 


HLZY-13 


0.02* 
0.22 
0.79 
0.16 
1.00 
0.14 
0.33 
0.72 
0.96 
0.24 
0.11 
0.01* 
0.00* 


HLZY-21 


0.36 
0.26 
0.02* 
0.21 
0.56 
0.00* 
0.33 
1.00 
0.61 
0.23 
0.02* 
0.09 
0.55 


Table 4 


Summary of results from analyses of molecular variance that indicate differences in 
genetic variation among and within distinct genetic groups of goldlined spinefoot (Sig- 
anus guttatus). Sequence data used in analyses are from specimens collected in the 
coastal waters of Vietnam during 2017—2020. The fixation index was 0.015. 


Sum of Meanof KEstimateof Percentage of 
Source of variation df squares squares variation variation (%) 


Among populations 2 23.756 11.878 0.060 1.5 
Within populations 344 705.852 2.051 0.114 98.5 
Total 693 1362.108 1.997 100.0 
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abundances and higher levels of biomass of target species 
were recorded inside MPAs as opposed to outside MPAs 
(Russ and Alcala, 2003; Maliao et al., 2009). However, the 
overall densities of goldlined spinefoot in MPAs were very 
low (with a mean of 0.51 fish/100 m7”), and these popula- 
tions may take a long time to recover to a high level of 
abundance. An effort to protect and conserve broodstocks 
at some MPAs with higher densities of goldlined spinefoot, 
especially at the Cu Lao Cham, Nha Trang Bay, and Phu 
Quoc MPAs, is an important measure that can be used by 
resource managers in the short term, and conservation 
at other MPAs can be useful in the longer term as well, 


because MPAs play an key role in enhancing local fish- 
eries outside MPAs through either net movement or net 
export of larvae from MPAs (Russ, 2002). 

The differentiation in population genetics of goldlined 
spinefoot between the central and southern groups is con- 
sistent with results of some previous studies, indicating 
that an oceanographic boundary limits the genetic disper- 
sal between the areas where these groups were found. 
Ablan et al. (2002) studied the genetic structure of 
threespot dascyllus and suggested that the South China 
Sea may be divided into 2 distinct parts on the basis of the 
differentiation in population genetics with geographical 
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Figure 5 


Visualization of the results from Bayesian analysis of the genetic 
structure of goldlined spinefoot (Siganus guttatus) collected from 
13 locations in the coastal waters of Vietnam during 2017-2020. 
Three distinct genetic groups were found in the sequence data from 
sampled fish. Bars represent the proportion of an individual fish’s 
genome that was assigned to the green, red, or blue genetic cluster. 
Specimens were collected from Con Co Island (CCO), Tam Giang— 
Cau Hai Lagoon (TAG), Hai Van—Son Cha (HVS), Thu Bon estuary 
(THB), the Cu Lao Cham Islands (CLC), Tra Bong estuary (TBO), 
the Ly Son Islands (LSO), Thi Nai Lagoon (TNA), Quy Nhon Bay 
(QNH), Nha Trang Bay (NTR), Nui Chua National Park (NCH), the 
Con Dao Islands (CDA), and the Phu Quoc Islands (PQU). 


boundaries located in south-central Vietnam close to Nha 
Trang Bay. A notable division in the genetic structure of 
goldlined spinefoot found in this study at the Phu Quoc 
Islands is in accordance with findings of the previously 
mentioned studies (Ablan et al., 2002; Dang et al., 2014). 
In addition, the results from recent studies based on anal- 
yses of communities of reef-building corals also indicate a 
separation in the distribution between the central and the 
southern groups, with a geographical boundary located at 
Dai Lanh Cape (at about 13°N) near the city of Tuy Hoa, 
and between the southern and southwestern groups at Ca 
Mau Cape (at about 8°N) at the southernmost tip of 
Vietnam (Vo, 2014; Vo and Nguyen, 2022). 

The discharge of the Mekong River (Dang et al., 2014), 
which acts as an oceanographic barrier to the dispersal of 
marine planktonic larvae, could be one explanation for the 
separation of reef fish communities between the south- 
western group in the Gulf of Thailand and the central 
and southern groups in the South China Sea. The recent 
results from oceanographic modeling indicate that a weak 
current flowed from the Gulf of Thailand (at the Phu Quoc 
Islands) to the Con Dao Islands, compared to a strong cur- 
rent from the Java Sea through the Karimata Strait, in 
Malaysia, up to the mouth of the Gulf of Thailand and over 
to the Con Dao Islands during the southwest monsoon 
(July-September), and that a weak current flowed from 
the Con Dao Islands to the Phu Quoc Islands, compared 
to a strong current from the Luzon Strait, between Tai- 
wan and the Philippines, east to Hainan Island, in China, 
and down along the coast of Vietnam through the Con 
Dao Islands, to Malaysia and the Karimata Strait during 
the northeast monsoon (January—March) (Li et al., 2019). 
These currents may limit the larval dispersal of goldlined 
spinefoot and other marine organisms. Our findings on 
the separation of the 3 genetic groups are similar to the 
findings of distinct populations of this species sampled in 
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Japan off Okinawa Island, Miyako Island, and 
Ishigaki Island on a geographical scale of about 
100—2000 km (Iwamoto et al., 2009). 

The spatial genetic structure found in our study 
for this species in the coastal waters of Vietnam 
could be related to differences in biological char- 
acteristics (Froukh and Kochzius, 2007). Larvae of 
goldlined spinefoot are in the pelagic stage, which 
lasts for about 24 d (Juario et al., 1985; Sheaves, 
1995; Nabhitabhata and Ikeda, 2014). This fac- 
tor may be the one that limits the dispersal of 
this species because larval dispersal of goldlined 
spinefoot is mainly driven by currents or tides 
(Kanashiro et al., 1999) and by seawater tempera- 
ture (Collins and Nelson, 1993). The numbers of 
alleles for groups of goldlined spinefoot were sim- 
ilar between sampling locations, with the excep- 
tion that the numbers for groups based on analysis 
of samples taken from specimens collected at the 
Phu Quoc Islands are twice as high as those for 
groups in other locations. This result indicates 
that there is distinguishable genetic diversity in 
this species at the Phu Quoc Islands. Therefore, 
the Phu Quoc Islands may be considered another functional 
location because of the differentiation in population genet- 
ics of this species. The differences in characteristics men- 
tioned for goldlined spinefoot at the Phu Quoc Islands may 
play important roles in the creation of the differentiation in 
population genetics among distinct groups. 


Conclusions 


The coastal lagoons and estuaries in Vietnam with seagrass 
beds (especially, Tam Giang—Cau Hai Lagoon, Nha Phu 
Lagoon, and Thu Bon estuary) supported high yield of juve- 
nile goldlined spinefoot, and the nearby MPAs (the Cu Lao 
Cham and Nha Trang Bay MPAs) had higher densities of 
adults of this species. There were 3 distinct genetic groups, 
with the central group including specimens from Con Co 
Island to Quy Nhon Bay, the southern group including 
individuals collected from Nha Trang Bay to the Con Dao 
Islands, and the southwestern group including fish from the 
Phu Quoc Islands in the Gulf of Thailand. Establishment 
of MPA networks on the basis of both ecological and social 
interconnectivity among crucial marine habitats, such as 
coral reefs, seagrass beds, and mangroves, is important 
because such efforts could help ensure that as many genetic 
units as possible are maintained for the overall resilience of 
fisheries. Creating networks of MPAs can further enhance 
the effectiveness of MPAs in protecting spawners in coral 
reefs and in increasing the survivorship of recruits in sea- 
grass beds and mangroves in coastal lagoons and estuaries. 


Resumen 


Las evaluaciones del estado de las especies objetivo de la 
pesca y sus vinculos ecolégicos con los habitats marinos 
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estan desempenando un papel importante en mejorar la 
efectividad de las areas marinas protegidas (AMP) y el 
manejo de otras zonas costeras. En este estudio, eval- 
uamos el rendimiento de juveniles y la abundancia del 
Ssigano rayas doradas (Siganus guttatus) utilizando el 
conocimiento local y censos visuales en aguas costeras 
de Vietnam. La estructura genética se examino medi- 
ante analisis de ADN de juveniles, colectados en pastos 
marinos y manglares en las principales lagunas y estu- 
arios, y de adultos colectados en arrecifes de coral. Se 
colectados mas de 38.6 millones de juveniles en 8 lagu- 
nas costeras y estuarios. Las lagunas costeras y los estu- 
arios con pastos marinos albergaban un alto volumen de 
juveniles, y las AMP cercanas presentaban mayores den- 
sidades de adultos. Se identificaron tres grupos genéti- 
cos distintos: un grupo en la parte central de la costa 
oriental, desde la isla de Con Co hasta la bahia de Quy 
Nhon; un segundo grupo en la costa sudoriental, desde 
la bahia de Nha Trang hasta las islas Con Dao, situadas 
en el extremo sur de Vietnam; y un tercer grupo en la 
costa sudoccidental, en las islas Phu Quoc, en el golfo de 
Tailandia. 
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Abstract—We examined the potential 
of North American river otters (Lontra 
canadensis) to buffer the expansion 
of the invasion by green crab (Carci- 
nus maenas) on the West Coast of the 
United States, documenting the diet of 
otters from scat remains on the Wa’atch 
and Tsoo-Yess Rivers, in Washington 
State, in 2018 and 2019. We tallied 
hard remains of prey and calculated fre- 
quency of occurrence, and we compared 
predation of the green crab to monthly 
values of catch per unit of effort for 
this crab species. North American river 
otters did not consume green crab in the 
Tsoo-Yess River and infrequently con- 
sumed green crab in the Wa’atch River 
(1.66% frequency of occurrence), likely 
because of the lower abundance of the 
green crab compared to the abundance 
of other prey in these rivers. Although 
our results indicate that North 
American river otters were not a biotic 
control of green crab, future studies on 
the population status of the green crab 
and North American river otter in both 
rivers and the long-term predator—prey 
dynamics could help to better gauge the 
potential for biotic resistance in popula- 
tions of the green crab. 
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The green crab (Carcinus maenas) 
is an invasive species known for its 
effects on eelgrass beds and juvenile 
bivalves (Cohen et al., 1995; Curtis 
et al., 2012; Howard et al., 2019). 
Green crab have been documented 
in the waters of Makah Bay on the 
Makah Indian Reservation in north- 
west Washington through removal 
trapping efforts, specifically efforts 
in the estuaries of the Wa’atch and 
Tsoo-Yess Rivers (Yamada?) (Fig. 1). 
The absence of green crab in nearby 
Neah Bay might be a result of compe- 
tition and predation by larger cancrid 
crab species, such as the red rock crab 


1 Yamada, S. B., C. Royer, S. Schooler, J. 
Fisher, A. Randall, C. Buffington, A. Stote, 
and A. Akmajian. 2022. Status of the 
European green crab, Carcinus maenas, 
in Oregon and Washington coastal estuar- 
ies: report for 2020 and 2021, 18 p. Report 
prepared for Aquatic Nuissance Species 
Project, Pacific State Marine Fisheries 
Commission. [Available from Aquat. Nuis- 
sance Species Proj., Pac. State Mar. Fish. 
Comm., 205 SE Spokane St., Ste. 100, 
Portland, OR 97202.] 


(Cancer productus) and adult Dunge- 
ness crab (Metacarcinus magister), 
both of which are more prevalent in 
Neah Bay than in the channels of the 
2 coastal rivers that enter Makah Bay 
(Jensen et al., 2007). On the western 
coast of Vancouver Island, in Canada, 
increasing abundances of Dungeness 
crab, red rock crab, or graceful crab 
(Metacarcinus gracilis) have resulted 
in reduced densities of green crab 
(Howard, 2019). Hence, biological 
resistance might be a principal fac- 
tor in mitigating the expansion of the 
invasion by green crab on the West 
Coast of the United States. 

The North American river otter 
(Lontra Canadensis) often forages in 
coastal environments in Washington 
State and is known to predate several 
marine crab species (Guertin et al., 
2010; Buzzell et al., 2014; Russell, 
2015). Little attention has been given 
to North American river otters 
and their potential to act as a biotic 
control for invasive species (Feltrop 
et al., 2016). Latrines (sites of marking 
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Figure 1 


Map of the sites (stars) where samples of scat from North American river otters (Lontra canaden- 
sis) were collected during the summer of 2018 in the Wa’atch River and during the spring and 
summer of 2019 in the Wa’atch and Tsoo-Yess Rivers. Collection sites were latrines, areas of 
marking and defecation by otters, located within the boundaries of the Makah Indian Reserva- 
tion in Washington State. Scat samples were used to examine the diet of North American river 
otters and the potential for predation by otters to control populations of the invasive green crab 
(Carcinus maenas). Map layer sources: WA State Parks GIS, ESRI, HERE, Garmin, SafeGraph, 
METI/NASA, USGS, EPA, USDA, NRCan, Parks Canada, GeoKye, Maxar, FAO, NOAA, Bureau 


of Land Management, NP. 


and defecation) of North American river otters have been 
observed along the lower Wa’atch and Tsoo-Yess Rivers 
during removal efforts targeting green crab (A. Akma- 
jian, personal observ.). Because of this apparent over- 
lap of habitats of North American river otters and green 
crab and the known predation of North American river 
otters on green crab in other regions (Mason and Mac- 
donald, 1980), we hypothesized that North American 
river otters are consuming green crab in both rivers and 
may act as a natural control for green crab in these 
estuaries. 

The first step in assessing the potential of North 
American river otters to serve as a biotic control of green 


crab was to determine 1) if North American river otters 
consume green crab and 2) if such consumption varies 
relative to time or abundance of green crab. We analyzed 
the hard remains found in scat from North American 
river otters to describe their diet and the relative impor- 
tance of prey items, focusing on consumption of green 
crab. We investigated differences in consumption of 
green crab between rivers, seasons (spring and summer 
of 2019), and years (summers of 2018 and 2019). Finally, 
to explore possible mitigation of the invasive green 
crab through their consumption by otters, we exam- 
ined the consumption of green crabs by North American 
river otters in relation to the availability of green crab 
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based on average catch rates during removal trapping 
efforts in 2019. 


Materials and methods 


The estuaries of the Wa’atch River and Tsoo-Yess River of 
the Makah Indian Reservation are located on the north- 
west coast of the Olympic Peninsula in Washington State 
(Fig. 1). The 2 rivers extend inland more than 16 km and 
up to 600 m in elevation, with the 4 river kilometers at 
the lowest elevations classified as brackish emergent 
tidal marsh (Heady et al., 2014). Active latrine sites used 
by North American river otters were identified along 
both the Wa’atch and Tsoo-Yess Rivers (Fig. 1). On the 
Wa’atch River, 1 latrine site was located on each of the 
north and south banks less than 1.5 km upstream from 
the river mouth. On the lower Tsoo-Yess River, 2 latrine 
sites were found on the east bank, approximately 0.8 and 
1.2 km upstream from the river mouth. The locations of 
the latrines in relation to the coast, bays, and one another 
are depicted in Figure 1. 

Prey identification was conducted by using analysis of 
undigested remains (bones and shells) found in scat sam- 
ples, a widely used method for documenting diet (Stenson 
et al., 1984; Day et al., 2015). Scat collection was completed 
every other week alongside green crab trapping efforts on 
the Wa’atch River in August and September 2018 and for 
a longer period in 2019 from April through September, but 
scat samples were collected on the Tsoo-Yess River only 
during May—September 2019. Samples of scat collected in 
the summer (July-September) from the latrines along the 
Wa’atch River were used to compare diet between years, 
but only samples collected in 2019 were used to compare 
diet between rivers and seasons (between spring, April— 
June, and summer). If latrines contained fewer than 15 
scat samples, collections were made on additional days in 
the following week. Protocols for the collection and clean- 
ing of scat were based on methods described for pinnipeds 
in Lance et al.” 

Prey items were identified to the lowest possible taxon 
but for brevity are reported herein only to the family level, 
apart from green crab and Dungeness crab. Color, claw, 
texture, and carapace morphology were all factors used to 
identify remains of crab. Fish remains were identified by 
using reference bones and otoliths housed at the NOAA 
Marine Mammal Laboratory in Seattle, Washington. We 
used frequency of occurrence (Trites and Joy, 2005) to 
examine and describe spatial and temporal variations in 
the diet of North American river otters, meaning differ- 
ences between the Tsoo-Yess and Wa’atch Rivers, between 
seasons (spring and summer), and between years. Mini- 
mum number of individuals was determined for crusta- 
ceans by using protocols adapted from Lance et al.” With 


2 Lance, M. M.,A. J. Orr, S. D. Riemer, M. J. Weise, and J. L. Laake. 
2001. Pinniped food habits and prey identification techniques 
protocol. Alsk. Fish. Sci. Cent. Proc. Rep. 2001-04, 29 p. [Avail- 
able from website, accessed November 2022. ] 


the exception of green crab, frequency of occurrence (FO) 
is reported for all prey groups with FO of at least 2% and 
was calculated as follows: 


where O, = 0 if taxon 7 is absent in fecal sample k and 1 if 
taxon i is present in fecal sample k; and 
s = total number of scat samples. 


We used catch per unit of effort (CPUE), defined as the 
number of green crab caught per trap set per day, as a proxy 
for abundance because no other information was available 
on populations of green crab in the study area. Trapping 
catch data from Makah Fisheries Management, Makah 
Tribe, were used in calculations of CPUE. Trap descrip- 
tions and protocols can be found in the summary report 
on trapping efforts in 2019 (Akmajian®). Average monthly 
values for CPUE of green crab and Dungeness crab were 
generated for each river to explore relative densities of 
prey. Although fish were also caught in traps, CPUE for 
other species in bycatch were not included in this study. 


Results and discussion 


In this study, 722 scat samples were collected in both riv- 
ers combined (Table 1). Green crab were present in sam- 
ples from the Wa’atch River (with FO of 1.66%) but absent 
in samples from the Tsoo-Yess River. Although green crab 
were scarce in the diet of North American river otters, 
green crab were observed in scat samples collected from 
both latrine sites along the Wa’atch River and on 7 sepa- 
rate dates, indicating that otters recognized green crab as 
prey (Suppl. Figure) (online only). Even so, green crab are 
elusive and aggressive, possibly making them more diffi- 
cult to capture and a less-favored prey type overall (Cohen 
et al., 1995). 

Fish composed the most frequently found and numerous 
group consumed by North American river otters in both 
rivers, a result that is consistent with those from other 
otter diet studies in coastal areas (Table 1) (Buzzell et al., 
2014; Russell, 2015). Among fish taxa, sculpins (Cottidae) 
made up the most frequently consumed prey group overall 
(65.37% FO), followed by gunnels (Pholidae) (57.20% FO) 
and righteye flounders (Pleuronectidae) (49.45% FO). 
Cancrid crabs (mostly Dungeness crab) composed the 
most consumed crustacean prey group (28.54% FO), but 
FO of cancrid crabs has varied in other studies (Table 1) 
(Guertin et al., 2010; Russell, 2015). 

The relative importance of prey groups in the diet of 
North American river otters was similar between 2018 
and 2019; however, most prey taxa and groups, including 
the green crab, had increased occurrence in scat samples 


3 Akmajian, A. 2020. European green crab trapping summary 
for the 2019 season, 30 p. [Available from Makah Fish. Manag., 
Makah Tribe, P.O. Box 115, Neah Bay, WA 98357.] 
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Table 1 


Frequency of occurrence (FO) for prey items found in samples of scat from North American river otters (Lontra 
canadensis) collected during the summer of 2018 in the Wa’atch River and during the spring and summer of 2019 
in the Wa’atch and Tsoo-Yess Rivers. Frequency of occurrence is reported for all prey groups that had FO >2%. The 
number of scat samples (7) is given for each collection period. 


Prey 


Teleostei 
Cottidae (sculpins) 
Pholidae (gunnels) 
Pleuronectidae (Righteye flounders) 
Embiotocidae (surfperches) 
Gasterosteidae (sticklebacks) 
Syngnathidae (pipefishes and seahorses) 
Salmonidae (salmons and trouts) 
Stichaeidae (pricklebacks) 
Teleostei (unidentified fish) 
Liparidae (snailfishes) 
Scorpaenidae (rockfishes) 
Decapoda 
Cancridae (rock crabs) 
Astacidae (crayfishes) 
Varunidae (shore crabs) 
Carcinus maenas (Portunidae) 
Miscellaneous 
Ranidae (true frogs) 


collected in 2019 (Table 1). With the exception of stickle- 
backs (Gasterosteidae) consumed by otters on the Tsoo-Yess 
River, the other 5 most frequently occurring prey groups 
increased in FO to some degree from spring to summer on 
both rivers. Seasonal change in the FO of green crab in the 
scat samples collected on the Wa’atch River was not appar- 
ent because of low consumption, where only 7 and 4 indi- 
viduals occurred in scat samples collected in spring and 
summer, respectively (Suppl. Figure) (online only). 
Dungeness crab had greater CPUE overall in trapping 
efforts in both rivers (1.40 individuals-set*-d-) than green 
crab (0.63 individuals-set™’-d). Such a difference coin- 
cides with the consumption of crabs by North American 
river otters, where both overall FO and the FO of Dunge- 
ness crab were higher than the FO of green crab (Table 1, 
Suppl. Figure) (online only). The CPUE of green crab was 
higher on the Tsoo-Yess River (0.84 individuals-set '-d~‘) 
than on the Wa’atch River (0.39 individuals-set *-d~4), but 
this difference was marginal compared to the greater over- 
all CPUE of Dungeness crab on both rivers during most of 
the trapping season. Although we used CPUE only as a 
proxy of prey abundance and availability, the difference 
between crab species in both CPUE and FO in the diet 
of North American river otters supports the hypothesis 


Wa’atch River Tsoo-Yess River 


2019 2019 


Spring Summer Spring Summer 
FO (%) FO (%) FO (%) FO (%) 
n=1385 w=179 n>125 n=144 


posited by others that this otter species consumes prey 
species in proportion to their relative abundance and 
availability (Ryder, 1955; Day et al., 2015). Therefore, low 
consumption of green crab was likely a result of the over- 
all lower abundance of green crab in both rivers compared 
with the abundance of other prey items. 

Several studies have documented the limitations of diet 
studies that use FO and minimum number of individuals 
(Crimmins et al., 2009; Tsukada et al., 2020). Scat collec- 
tion is subjective and can lead to pseudoreplication, where 
there is a strong likelihood that multiple scat samples from 
a single collection are from the same individual, leading 
to an overestimate of frequently consumed species and an 
underestimate of infrequently consumed species (Tsukada 
et al., 2020). Estimating species abundance by using 
CPUE has limitations because of the high variability in 
trapping efficacy; in future studies, an aim should be to 
use methods, such as mark-recapture techniques, that are 
more accurate (Munch-Petersen et al., 1982, Bergshoeff 
et al., 2018; Bernier et al., 2020). Although it is important 
to consider these biases, the objective of this study was to 
detect green crab in the diet of North American river otters 
and to describe important prey types, which was possible 
because of the timeframe of our study. 
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On the basis of our results, it is unlikely that North 
American river otters currently serve as a natural control 
of green crab populations that are at low densities. How- 
ever, in recent years, population densities of green crab 
have continued to increase in Makah Bay and elsewhere 
along the coastline of Washington. Predation of green 
crab by North American river otters will likely increase 
as green crab become more abundant and available, but 
gauging potential effects will require estimating popula- 
tion densities for both North American river otters and 
green crab. As such, future studies should focus on under- 
standing predator—prey interactions involving green crab 
and other likely predators, such as adult Dungeness crab, 
that might cumulatively buffer a higher density of green 
crab and further mitigate effects of green crab invasions 
on sensitive coastal habitats along the West Coast. 


Resumen 


Examinamos el potencial de las nutrias de rio norteamer- 
icanas (Lontra canadensis) para amortiguar la expansi0n 
de la invasién del cangrejo verde (Carcinus maenas) en 
la costa oeste de Estados Unidos, documentando la dieta 
de las nutrias a partir de restos de excrementos en los 
rios Wa’atch y Tsoo-Yess, en el estado de Washington, en 
2018 y 2019. Contamos los restos duros de presas, calcu- 
lamos la frecuencia porcentual de ocurrencia, y compar- 
amos la depredacién del cangrejo verde con los valores 
mensuales de captura por unidad de esfuerzo para esta 
especie de cangrejo. Las nutrias de rio norteamericanas 
no consumieron cangrejos verdes en el rio Tsoo-Yess y 
consumieron con poca frecuencia cangrejos verdes en el 
rio Wa’atch (1.66% de frecuencia de ocurrencia), probable- 
mente debido a la menor abundancia del cangrejo verde 
en comparacioén con la abundancia de otras presas en 
dichos rios. No obstante que nuestros resultados indican 
que las nutrias de rio norteamericanas no constituyen 
un control biologico del cangrejo verde, futuros estudios 
sobre el estado de las poblaciones del cangrejo verde y 
nutria de rio norteamericana en ambos rios y la dinaémica 
depredador-presa a largo plazo podrian ayudar a calibrar 
mejor el potencial de resistencia bidética de las poblaciones 
del cangrejo verde. 
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Abstract—The results of this study 
provide additional evidence of ovar- 
ian masculinization in the Atlantic 
croaker (Micropogonias undulatus) in 
and around the hypoxic zone in the 
Gulf of Mexico and the first evidence 
of ovarian masculinization in the spot 
(Leiostomus xanthurus) and _ bigeye 
searobin (Prionotus longispinosus). In 
some specimens of Atlantic croaker and 
spot, we observed this phenomenon to 
a degree that has not been reported in 
either species. We also found evidence 
of smaller gonads in spot from hypoxic 
sites compared with the gonads in spot 
from sites with normal oxygen concen- 
trations. These findings indicate that 
Ovarian masculinization and _ repro- 
ductive impairment associated with 
hypoxia exposure affect a wider range 
of species than was previously known, 
and they highlight the importance of 
continued research on this topic. 
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In the last century, hypoxia has emerged 
as a significant threat to coastal eco- 
systems. Hypoxia, or hypoxic condi- 
tions, refers to low levels of dissolved 
oxygen in bodies of water. Coastal 
hypoxia, which was first observed in the 
1930s (Diaz and Rosenberg, 2008), has 
increased 3-fold during the last 30 years 
(Diaz and Rosenberg, 2008; Vaquer- 
Sunyer and Duarte, 2008). The range of 
dissolved oxygen values used to define 
an area as hypoxic is variable, and the 
term hypoxia conventionally indicates 
dissolved oxygen values below 2.0 mg/L 
(Rabalais et al., 2002; Vaquer-Sunyer 
and Duarte, 2008). Although hypoxic 
events are known to occur naturally 
throughout the world (Helly and Levin, 
2004), the global increase in coastal 
hypoxia can mostly be attributed to the 
increased use of synthetic nitrogen fer- 
tilizers since their creation in the 1940s 
(Diaz and Rosenberg, 2008; Galloway 
et al., 2008). Increased nitrogen runoff 
from fertilizers can result in large algal 
blooms in coastal waters. When algal 
cells die, their decomposition consumes 
dissolved oxygen in the water column. 
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This process can result in hypoxia if the 
water column is stratified, given that 
subsurface waters are unable to mix 
with oxygen-rich surface waters because 
of differences in water density. 

Hypoxia has wide-ranging ecological 
and physiological effects on marine life. 
Exposure to hypoxic conditions alters 
the behavior (Eriksson and Baden, 
1997) and interspecific interactions of 
marine organisms (Sandberg, 1997). 
Hypoxic conditions have also been asso- 
ciated with changes in community com- 
position (Dauer, 1993; Pihl, 1994) and 
biomass (Rosenberg and Loo, 1988). 
Community recovery following epi- 
sodes of hypoxia is highly variable and 
dependent upon a variety of environ- 
mental and biotic factors (Wu, 2002), 
including community complexity and 
the severity and duration of hypoxic 
conditions (Boesch and _ Rosenberg, 
1981). At the organismal level, endo- 
crine disruption (Wu et al., 2003) and 
reproductive impairment (Zhou, 2001) 
are 2 of the most well-documented 
physiological effects of hypoxia expo- 
sure in fish. 
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Fish exposed to hypoxic conditions tend to have lower 
values for the gonadosomatic index (GSI), an indicator of 
reproductive condition, than fish exposed to normal oxy- 
gen concentrations, as has been observed in the Gulf killi- 
fish (Fundulus grandis) (Cheek et al., 2009), common carp 
(Cyprinus carpio) (Zhou, 2001), and Atlantic croaker (Micro- 
pogonias undulatus) (Wu et al., 2003; Thomas et al., 2007; 
Thomas and Rahman, 2012). Both episodic (diel) and chronic 
hypoxia have been shown to reduce the reproductive capac- 
ity of estuarine fish species (Landry et al., 2007; Cheek et al., 
2009). Chronic exposure to hypoxic conditions has been fur- 
ther linked to masculinization in female fish (Landry et al., 
2007; Cheek et al., 2009; Thomas and Rahman, 2012) and 
to reduced sperm production and viability in male Atlantic 
croaker (Thomas et al., 2015). Masculinization of female 
Atlantic croaker resulting from chronic exposure to hypoxia 
has been documented in both field and laboratory studies 
(Thomas et al., 2007; Thomas and Rahman, 2012). 

In this study, we assessed whether the well-established 
reproductive effects of seasonal hypoxia on Atlantic croaker 
in the Gulf of Mexico (GOM) (Thomas and Rahman, 2012) 
extend to other cohabiting species of fish, specifically the 
spot (Leiostomus xanthurus) and bigeye searobin (Prion- 
otus longispinosus). Both the Atlantic croaker and spot 
are members of the family Sciaenidae (order Perciformes). 
The bigeye searobin belongs to the family Triglidae (order 
Scorpaeniformes). Finding evidence of reproductive abnor- 
malities at the family and order levels would demonstrate 
that the effects of seasonal hypoxia in the GOM are more 
widespread taxonomically, with broader implications for 
fisheries management, than previously thought. 


Materials and methods 
Sample collection and study area 


Fish were collected aboard the NOAA Ship Oregon IT during 
legs I and IT (June-August) of the Southeast Area Monitor- 
ing and Assessment Program (SEAMAP) summer ground- 
fish surveys in 2016, 2017, and 2019 (Fig. 1). Sampling sites 
were scattered throughout the area in which hypoxia season- 
ally develops in the GOM, ranging from the coast of Texas 
to waters off Louisiana. Samples were collected as described 
in the SEAMAP sampling protocol (GSMFC°) at survey sites 
with depths <30 m because hypoxia in the GOM occurs less 
frequently at depths below 30 m. Whole specimens caught 
during trawl surveys were fixed in a 10% formalin solution 
and transferred to ethanol for storage. Dissolved oxygen val- 
ues were recorded immediately prior to trawling, according 
to SEAMAP protocol for collecting data on water conditions, 
by using a shipboard array equipped with dissolved oxygen, 
conductivity, temperature, and depth sensors and water bot- 
tles for measuring other water properties (Table 1). A site was 


1 GSMFC (Gulf States Marine Fisheries Commission). 2019. 
SEAMAP operations manual for trawl and plankton surveys, 
63 p. Gulf States Mar. Fish. Comm., Ocean Springs, MS. [Avail- 
able from website.] 


defined as hypoxic if dissolved oxygen was 2.0 mg/L or lower, 
a level that is consistent with the recent literature on hypoxia 
(Vaquer-Sunyer and Duarte, 2008). 

Female specimens of Atlantic croaker, spot, and bigeye 
searobin were reproductively immature at the time of sam- 
pling, as they did not contain any oocytes past the previ- 
tellogenic stage (Anderson et al., 2018). Male specimens of 
Atlantic croaker and spot were also reproductively imma- 
ture, as no spermatids or spermatozoa were present. Some 
male bigeye searobin were reproductively mature. The nat- 
ural history of Atlantic croaker and spot and histological 
examination of gonads of bigeye searobin strongly indicate 
that specimens would have reproductively matured within 
the year, making any evidence of reproductive abnormal- 
ities indicative of reproductive impairment. It is import- 
ant to note that, although the majority of specimens were 
reproductively immature, the sampled fish were not juve- 
niles and the term immature is meant to describe the state 
of the gametes within the gonads. 

We obtained reference samples of Atlantic croaker and 
spot from the Royal D. Suttkus Fish Collection at the Tulane 
University Biodiversity Research Institute (TUBRD. Refer- 
ence samples consisted of whole, formalin-fixed specimens 
from which gonads were dissected. We selected reference 
samples collected during the summer (June—September) at 
either the Chandeleur Islands or the Rigolets, both located 
in Louisiana. We selected samples collected between 1955 
and 1969 because we wanted specimens from a time before 
the 1980s, when eutrophication dramatically increased in 
the western GOM, as evidenced by increased carbon depo- 
sition (Rabalais et al., 2007). The Chandeleur Islands and 
the Rigolets are not known to be hypoxic during the sum- 
mer; therefore, both the location and time of the collection 
of our references samples makes it unlikely that reference 
fish were exposed to hypoxic conditions prior to sampling. 


Sample processing 


Subsamples were selected from male and female Atlan- 
tic croaker, spot, and bigeye searobin collected from both 
hypoxic sites and sites with normal oxygen levels. These 
selected specimens were wet weighed (to the nearest 
0.001 g) and measured for total length (to the nearest 
millimeter), and then they had their gonads removed and 
wet weighed. The gonads were mounted, sectioned, and 
stained with hematoxylin and eosin. Only a subsample of 
these gonads were histologically processed and examined. 


Gonadosomatic index 


The GSI was calculated by using this formula: 


Weight 
§ gonads x 100, 
(Weights, — Weight gonads) 


= the wet weight of both gonads post 
fixation; and 

Weight, = the wet weight of the whole fish post 

fixation. 


where Weight 


gonads 
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Figure 1 


A map of the sites in the northern Gulf of Mexico where Atlantic croaker (Micropogonias undulatus), spot 
(Leiostomus xanthurus), and bigeye searobin (Prionotus longispinosus) were sampled in 2016, 2017, and 
2019. All sampling for this study was conducted during Southeast Area Monitoring and Assessment Program 
(SEAMAP) summer groundfish surveys aboard the NOAA Ship Oregon IT at sites with depths <30 m scattered 
throughout the hypoxic zone that develops seasonally in the Gulf of Mexico. Reference samples, collected 
during summer between 1955 and 1969 at sites known not to have hypoxic conditions, were obtained from the 
Royal D. Suttkus Fish Collection of the Tulane University Biodiversity Research Institute. Numerals next to 
the symbols for sites indicate the station numbers assigned by SEAMAP. Note that the shoreward boundaries 
of the polygons for hypoxic areas reflect the boundaries of the hydrographic survey used to generate the poly- 
gons and do not necessarily represent the extent of the hypoxic zone in a given year. 


In order to compare the GSI of fish of varying sizes, a 
residual index of the GSI values was computed. The resid- 
ual index has been shown to be an effective way to reduce 
the influence of body size in analyses of organismal condi- 
tion (Jakob et al., 1996) and has been proven to be effective 
in the analysis of fish morphometric data (Reist, 1985; 
Jakob et al., 1996). The residual index is the difference 
between the observed GSI value and the predicted GSI 
value, which is determined by using a linear trendline gen- 
erated with pooled data (for specimens from both hypoxic 
and normal oxygen concentrations). Positive values for the 
residual index indicate that the observed GSI is greater 
than the predicted GSI, whereas negative values indicate 
that the observed GSI is less than the predicted GSI. 


Overall body condition 


Overall body condition was calculated as wet weight 
divided by total length. These condition index values for 


specimens were compared between hypoxic and normal 
sites by using the residual index (polynomial trendline) 
of the length—weight relationship based on data pooled for 
males and females, under the assumption that there is no 
sexual dimorphism in any of the species included in this 
study. In order to compare both sexes simultaneously, the 
weight used for the overall condition analysis was the wet 
weight of each fish minus the weight of its gonads. 


Statistical treatment 


Mean residual GSI and condition index values were com- 
pared between fish from hypoxic sites and those from 
sites with normal oxygen levels by using sites as repli- 
cates. The residual indices of both overall body condition 
and GSI for both male (hypoxic: n=31; normal: n=28) and 
female (hypoxic: n=36; normal: n=31) spot and female 
Atlantic croaker (hypoxic: n=37; normal: n=44) were nor- 
mally distributed (Shapiro—Wilk test: P>0.05) (Shapiro 
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Table 1 


Recorded dissolved oxygen (DO) values and depths at the 
sites sampled for this study in the northern Gulf of Mexico 
during 2016, 2017, and 2019. The site numbers correspond 
to the numbers assigned to “Pascagoula” stations by the 
Southeast Area Monitoring and Assessment Program. 


DO Depth 
(mg/L) (m) 


Site no. Year 


15 2016 5.7 22.4 
32 2016 1 15.4 
98 2016 1.4 19.6 
99 2016 iE 18.7 
108 2016 2.3 22.8 
ir 2016 0.0 15.0 
118 2016 0.9 14.4 
119 2016 1.0 15.8 
120 2016 5.1 15-9 
177 2016 4.1 26.4 
179 2016 2.2 23.8 
198 2016 2.0 18.1 
199 2016 1.8 18.8 
28 2017 5.8 23.6 
29 2017 6.3 19.4 
36 2017 5.4 18.0 
37 2017 2.5 21.6 
=o 2017 1.4 20.1 
45 2017 3.4 18.8 
47 2017 3.2 19.3 
48 2017 5.9 15.6 
Fak 2017 1.4 17.8 
46 2019 3.0 33.0 
63 2019 0.5 14.0 
64 2019 0.5 13.0 
72 2019 1.5 10.0 
98 2019 1.5 26.0 


and Wilk, 1965) and had statistically indistinguishable 
variances (Bartlett’s test: P>0.05) (Bartlett, 1937). There- 
fore, we used a 2-way analysis of variance (ANOVA), and 
a significance level (a) of 0.05, to compare the residual 
indices of these samples. The residual GSI for male 
Atlantic croaker was not normally distributed (Shapiro— 
Wilk test: P<0.05); therefore, we used a Kruskal—Wallis 
test, with a o of 0.05, for these data (hypoxic: n=48; nor- 
mal: n=26) (Kruskal and Wallis, 1952). 


Results 
Ovarian masculinization 


Evidence of masculinization was found in 6 of the 17 ova- 
ries from Atlantic croaker that were examined in this 
study, with masculinized ovaries observed only for speci- 
mens collected at hypoxic sites. In the masculinized ova- 
ries from Atlantic croaker, testicular tissue was found 
primarily in connective tissue along the ovary periph- 
ery (Fig. 2). The testicular tissue consisted of clusters 


of spermatids within connective tissue, similar to the 
tissue conditions used in determinations of masculin- 
ization made in previous studies (Thomas and Rahman, 
2012) (Fig. 2). The spermatogenic tissue observed within 
ovarian tissue resembled spermatogenic cysts of lobular 
type testes. One specimen, the gonads of which exter- 
nally resembled testes, had gonadal tissue that is best 
described as transitional (i.e., an ovary that was transi- 
tioning to a testis). The gonads of this specimen organiza- 
tionally resembled testes; however, degenerated oocytes 
were present and surrounded duct-like structures. Sper- 
matogenic tissue was not present. The ovaries of refer- 
ence specimens contained no evidence of spermatogenic 
tissue or ovarian anomalies. 

The masculinized ovaries from spot (Fig. 3) had a mas- 
culinization profile similar to that of the masculinized ova- 
ries from Atlantic croaker. Masculinization in peripheral 
connective tissue was found in 3 of the 7 ovaries from spot 
that were examined in this study, with one of the ovaries 
coming from a specimen collected at a site with normal 
oxygen concentrations. 

Two of the spot we examined contained ovaries that 
are best described as in a transition from female to male 
tissue (Fig. 3). These transitional gonads had structures 
more like a typical testis, with ducts surrounded by con- 
nective tissue that was interspersed with spermatogenic 
tissue. One of these transitional gonads contained numer- 
ous atretic follicles and a large lumen filled with what 
looked like ovum remnants and large melano-macrophage 
centers (Agius and Roberts, 2003). The other transitional 
gonad contained atretic follicles and a testicular-like struc- 
ture but lacked ovum remnants and melano-macrophage 
centers. The ovaries of the reference specimens used in 
this study appeared more organized compared with the 
gonads of specimens collected in 2016, and no evidence of 
spermatogenic tissue was found within them. 

Masculinization in bigeye searobin manifested as sper- 
matogenic cells in various stages of development inter- 
spersed among oocytes and peripheral tissue (Fig. 4). 
Evidence of this phenomenon was found in 1 of the 7 ova- 
ries from bigeye searobin that were examined in this study. 


Gonadosomatic index and body condition 


The mean residual index of overall body condition based 
on pooled data for female and male Atlantic croaker 
was significantly higher for fish from hypoxic sites than 
for fish from sites with normal oxygen concentrations 
(Kruskal—Wallis test: P<0.05) (Figs. 5 and 6.). The resid- 
ual GSI of both male and female Atlantic croaker from 
hypoxic sites did not significantly differ from that of 
males and females from sites with normal levels of oxy- 
gen (males: P=0.65, df=7, f=0.21; females: P=0.77, df=11, 
f=0.08) (Figs. 5 and 6). 

The mean residual index of overall body condition for 
female and male spot combined did not significantly dif- 
fer between specimens from hypoxic sites and specimens 
from normal sites (P=0.12, df=14, f=0.60) (Figs. 7 and 8). 
The mean residual GSI of female spot from hypoxic sites 
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Figure 2 


Microphotographs of the gonads of Atlantic croaker (Micropogonias undulatus) collected from the 
northern Gulf of Mexico, showing the effects of hypoxia on gonadal development. (A) This ovary 
of a reference specimen, acquired from the Tulane University Biodiversity Research Institute and 
sampled between 1955 and 1959, has previtellogenic oocytes (Oo) and no signs of spermatogenic 
tissue or oocyte degradation. (B) In this mature testis of a male, from a control treatment of a 
laboratory study conducted in 2020 (senior author, in prep.), spermatocytes (Sc) in various stages 
of development are visible, as are fully developed spermatids (S). No cells in this mature testis 
resemble oocytes. (C) In this section of an ovary, from a specimen collected in 2016 at a site ina 
hypoxic area (site 99, which had a dissolved oxygen [DO] level of 1.1 mg/L), both spermatocytes are 
in various stages of development. (D) This transitional gonad, from a specimen externally classi- 
fied as male and sampled in a hypoxic zone (in an area with a DO level of 1.4 mg/L) in 2016, has 
clear ducts (D) surrounded by degrading or degenerating oocytes. The structure resembles a testis, 
with reproductive cells surrounding duct-like structures. Degenerating oocytes in this gonad are 
larger than and visually distinct from spermatocytes in the mature testes shown in panel C. 


was significantly lower (ANOVA: P<0.05) than that of 
females from sites with normal oxygen concentrations 
(Figs. 7 and 8), and male spot from hypoxic sites also had 
significantly lower GSI values than males from normal 
sites (ANOVA: P<0.05). 


Discussion 

Scope and possible causes of ovarian masculinization 

Prior to this study, the only species from the GOM for 
which evidence of ovarian masculinization had been 


reported was the Atlantic croaker (Thomas and Rahman, 
2012). This study is the first to report this phenomenon in 


2 additional GOM species, the spot and bigeye searobin, 
with 2 of 7 subsampled ovaries from spot (Fig. 3) and 1 of 
7 subsampled ovaries from bigeye searobin (Fig. 4) con- 
taining male reproductive cells. Additionally, such evi- 
dence was found in 6 of 17 ovaries of Atlantic croaker 
examined in this study (Fig. 2). 

The cause of the ovarian masculinization found in this 
study cannot be determined. However, hypoxia is a likely 
possibility. It is well known that estrogen, the synthesis of 
which is regulated by the enzyme aromatase, is a primary 
factor in the sex determination of fish. Guiguen et al. 
(2010) discussed results of a number of experiments that 
indicate that suppression of estrogen synthesis (aromatase 
inhibition) in fish can result in the formation of testicular 
tissue. Aromatase inhibitors have been shown to induce 


Cyrana et al.: Ovarian masculinization and reproductive impairment in 3 species of groundfish 


Figure 3 


Microphotographs of the gonads of spot (Leiostomus xanthurus) collected from the northern Gulf 
of Mexico, showing the effects of hypoxia on gonadal development. (A) This ovary of a reference 
specimen, acquired from the Tulane University Biodiversity Research Institute and sampled in 
1969, has previtellogenic oocytes (Oo) and no evidence of spermatogenic tissue or oocyte degra- 
dation. (B) In this ovary of a specimen collected in 2016, at a site adjacent to an area of hypoxia 
(site 177, which had a level of dissolved oxygen [DO] of 4.1 mg/L), oocytes are visible, as are clus- 
tered spermatozoa (Sz). (C) This transitional gonad, from a specimen sampled in 2016 at a site 
in a hypoxic zone (site 198, which had a DO level of 2.0 mg/L), has testicular features, like ducts 
(D) and a preponderance of connective tissue, and evidence of previous oocyte degradation, such 
as degraded oocyte material (DOM) and melano-macrophage aggregations (Mm). (D) In the same 
specimen featured in panel C, shown here at a higher magnification, degraded or degenerating 


oocytes are embedded in the connective tissue. Testicular-like ducts are also visible. 


masculinization in several species of fish, including the 
naturally hermaphroditic honeycomb grouper (Epi- 
nephelus merra) (Bhandari et al., 2004), red grouper 
(E. morio) (Guiguen et al., 2010), and European seabass 
(Dicentrarchus labrax) (Navarro-Martin et al., 2009) and 
in the naturally gonochoristic common carp (Ogawa et al., 
2008). Long-term aromatase suppression has even been 
reported to convert functional ovaries into functional tes- 
tes in tilapia (Sun et al., 2014). Hypoxic conditions have 
been shown to function as an aromatase inhibitor in 
Atlantic croaker, decreasing mRNA expression and aro- 
matase activity in fish exposed to hypoxia both in the wild 
and in the lab (Thomas and Rahman, 2012). 

It is likely that the ovarian masculinization observed 
in fish from hypoxic sites in our study could be a conse- 
quence of hypoxia-induced inhibition of aromatase and the 


concomitant suppression of estrogen synthesis. The well- 
documented connection between aromatase suppression 
and the development of intersex gonads in fish, coupled 
with the broad geographic extent of hypoxia encountered 
during the sampling for this study, strongly indicates that 
the ovarian masculinization observed in this study was 
hypoxia related. 

The presence of intersex specimens from sites with nor- 
mal oxygen levels is a novel discovery; however, the rela- 
tively low sample size prevents us from offering a definitive 
interpretation. Specifically, the normal sites that yielded 
intersex specimens were immediately adjacent to areas of 
hypoxia (Fig. 1) and differed by only 0.2 mg/L in the level 
of dissolved oxygen, with the concentration at these nor- 
mal sites falling right on the line between our definition of 
hypoxic and normal conditions. Although this closeness in 
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Figure 4 


Microphotographs of ovaries of bigeye searobin (Prionotus 
longispinosus) collected in the northern Gulf of Mexico in 
2016. (A) This ovary has no signs of masculinization or 
oocyte degradation: previtellogenic oocytes (Oo) are visible 
throughout the ovary. (B) In this ovary, there is evidence of 
masculinization: both oocytes and spermatocytes (Sc) are 
in various stages of development, and spermatozoa (Sz) 
are present. 


location and dissolved oxygen could indicate either that 
the intersex specimens had passed through the hypoxic 
area and were subsequently masculinized or that the dis- 
solved oxygen concentration at the normal site was low 
enough to induce the development of intersex gonads, 
there are not enough data to identify a causal relation- 
ship. Sampling numbers were not sufficient to explain 
the source of ovarian masculinization in fish from normal 
sites; therefore, further sampling efforts and laboratory 
exposure experiments are needed to fully explain this phe- 
nomenon. Similarly, it is unclear what effect the masculin- 
ization observed in this study would have had on fish if 
they had reached reproductive maturity. Although the 
molecular pathway involved in ovarian masculinization in 
fish is well understood, the parameters that trigger this 


pathway and how they affect fish at different stages of 
reproductive development is not well understood and 
should be further studied. 

Our findings highlight the importance of continued 
work in this area, particularly in the development of inter- 
sex gonads in additional fish species and the exposure 
parameters (e.g., duration of exposure, dissolved oxygen 
level, and reproductive stage at time of exposure). A bet- 
ter understanding of the exposure parameters that trigger 
the development of intersexuality in fish would help not 
only elucidate the findings of this study but also better 
inform management decisions regarding resources in the 
GOM and in other areas affected by hypoxia. 

Our observation of ovarian masculinization in Atlantic 
croaker and 2 other groundfish species in the GOM 
extends the taxonomic breadth of this phenomenon, which 
was first reported by Thomas and Rahman (Thomas 
et al., 2007). If our finding of testicular tissue in the ova- 
ries of wild-caught specimens is not related to hypoxia, 
the phenomenon either is due to some other environ- 
mental factor or it represents naturally occurring back- 
ground intersexuality. Bahamonde et al. (2013) reviewed 
this phenomenon in great detail, with his work includ- 
ing studies of gonochoristic species in which intersex 
individuals were found to occur naturally (Bahamonde 
et al., 2013). In one such study involving the freshwa- 
ter species smallmouth bass (Micropterus dolomieu), 
significantly more occurrences of intersex fish during 
the prespawning period were observed than during the 
reproductive period (Blazer et al., 2007). The spawning 
period of Atlantic croaker in the GOM begins as early as 
October and ends as late as June (Parker, 1971), a pro- 
tracted spawning period. Results from research on bigeye 
searobin in the GOM indicate that this species also has a 
protracted spawning period, from January through June 
(Hoff, 1992). The spot has been reported to be a winter 
spawner in the GOM, with the exact timing of spawning 
contingent on longitude (Parker, 1971). Even though our 
ovarian samples were taken from fish collected primarily 
in the prespawning period for their species, it is unlikely 
that the evidence of masculinization we observed was the 
result of the inherent intersex characteristics described 
by Bahamonde et al. (2013). 

The vast majority of examples of naturally occurring 
intersex individuals in gonochoristic fishes outlined by 
Bahamonde et al. (2013) are male fish with testicular tis- 
sue that contained oocytes, given that the study was 
focused primarily on the process of male feminization. 
The feminization of testicular tissue indicates aromati- 
zation of androgens into estrogens—the opposite hor- 
monal mechanism from that presumed to occur in ovarian 
masculinization, making it incompatible with our inter- 
pretation. Additionally, our observation of heavily mascu- 
linized specimens, including a specimen of spot best 
characterized as in the process of undergoing a sex 
change, with degenerating oocytes, indicates that speci- 
mens were transitioning from female to male, similar to 
the transition of ovaries to functional testes in tilapia as 
a result of aromatase suppression (Sun et al., 2014). 
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(A) Weight-—length relationship for Atlantic croaker (Micropogonias undulatus), males and females combined, and the 
relationship between length and the gonadosomatic index for (B) female and (C) male Atlantic croaker collected in the 
northern Gulf of Mexico during 2016, 2017, and 2019. All specimens were caught off the coasts of Texas and Louisiana 
at sites in hypoxic areas and sites with normal levels of oxygen. The solid trend lines were calculated by using data 
pooled for fish from both hypoxic and normal sites, and the dashed and dotted trend lines represent data only for fish 
from hypoxic sites or only for fish from normal sites, respectively. Also provided are the coefficient of determination 


(r?) in each graph and the equation for the solid trend line in each of the bottom graphs. 


Therefore, although the occurrence of naturally intersex 
individuals in otherwise gonochoristic fishes should be 
seriously considered, the low number of examples of nat- 
urally masculinized individuals in the literature, com- 
pared to the high degree of masculinization in some of 
our specimens, as well as the lack of any observable 
intersexuality in our reference samples indicate that nat- 
urally occurring intersexuality is not a likely explanation 
for our observations. 


Manifestation of masculinization 


In previous reports of masculinization in Atlantic croaker, 
the process is described as “spermatogenic cyst-like struc- 
tures” (Thomas and Rahman, 2012). This description fits 
the majority of masculinized ovaries from Atlantic croaker 
and the masculinized specimens of spot from a site in this 
study that had normal oxygen concentrations but was 
adjacent to the area with hypoxic conditions. In addition to 
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Mean residual gonadosomatic index (GSI), mean residual overall body condition index, and box plots of GSI values 
for Atlantic croaker (Micropogonias undulatus) caught at hypoxic sites and sites with normal oxygen concentrations 
in the northern Gulf of Mexico in 2016, 2017, and 2019. (A) Mean residual GSI and mean residual index of overall 
body condition, calculated as wet weight divided by total length, for females and males (analysis of variance: P<0.05; 
hypoxic sites: n=9; normal sites: n=11) are provided. Sampling sites served as replicates for all analyses. Error bars 
indicate standard errors of the means. For (B) female and (C) male Atlantic croaker, box plots present residual GSI 
values by site and site type. Numerals along the x-axis indicate the numbers assigned to sampling sites by the South- 
east Area Monitoring and Assessment Program and the year during which specimens were collected. The upper and 
lower parts of each box represent the first and third quartiles (the 25th and 75th percentiles), and the horizontal line 
is the median. The whiskers extend above and below the box no more than 1.5 times the interquartile range. Data 
points inside and outside the boxes represent residual GSI values for individual fish. 


observing this more typical manifestation of masculiniza- 
tion in sampled fish, we also examined 3 specimens that 
were undergoing this phenomenon to a degree that has yet 
to be described in wild specimens from the GOM. 

In 2 specimens of spot and 1 specimen of Atlantic 
croaker, we observed masculinization to the extent that 
their gonads were in transition from female to male tissue, 
as these specimens had degenerating or regressing female 
reproductive cells and a testes-like overall structure. 
The gonads of these specimens contained well-defined 
ducts that were surrounded by highly degraded oocytes. 


The gonad of one of the spot specimens also contained 
a large duct or lumen space filled with material consis- 
tent with follicular atresia and melano-macrophages. The 
other specimen of spot had similarly degraded oocytes 
and well-defined male tissue; however, the gonad of this 
individual lacked a large lumen and melano-macrophage 
aggregations. 

The evidence of oocyte degradation combined with the 
well-defined testicular structure and spermatogenic tis- 
sue in these specimens indicate masculinization beyond 
what has been observed previously in fish from the GOM. 
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Figure 7 


(A) Weight-—length relationship for spot (Leiostomus xanthurus), males and females combined, and the relationship 
between length and the gonadosomatic index for (B) female and (C) male spot collected in the northern Gulf of Mexico 
during 2016, 2017, and 2019. All specimens were caught off the coasts of Texas and Louisiana at sites in hypoxic areas 
and sites with normal levels of oxygen. The solid trend lines were calculated by using data pooled for fish from both 
hypoxic and normal sites, and the dashed and dotted trend lines represent data only for fish from hypoxic sites or only 
for fish from normal sites, respectively. Also provided are the coefficient of determination (r”) in each graph and the 


equation for the solid trend line in each of the bottom graphs. 


The characterization of these individuals as females 
undergoing masculinization, as opposed to males undergo- 
ing feminization, is supported by work on tilapia by Sun 
et al. (2014) in which aromatase suppression led to a com- 
plete sex change in females, a change characterized by a 
degradation of oocytes and replacement of the oocytes by 
spermatogenic cells (Sun et al., 2014). Histological obser- 
vations in naturally hermaphroditic fish transitioning 
from female to male provide similar evidence of oocyte 


breakdown prior to the transition from female to male 
(Wang et al., 2017; Maxfield and Cole, 2019). Centers of 
melano-macrophages have been considered pathological 
signs of environmental degradation and hypoxia in fish 
(Agius and Roberts, 2003), including individuals of species 
such as the Atlantic croaker in the GOM (Fournie et al., 
2001), albeit the aggregations were found in the spleen. 
The more structurally complex nature of the examples of 
masculinization we observed in this study, coupled with 
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Mean residual gonadosomatic index (GSI), mean residual overall body condition index, and box plots for spot (Leiosto- 
mus xanthurus) caught at hypoxic sites and sites with normal oxygen concentrations in the northern Gulf of Mexico 
in 2016, 2017, and 2019. (A) Mean residual GSI for females (ANOVA: P<0.05; hypoxic sites: n=6; normal sites: n=8) 
and males (ANOVA: P< 0.05; hypoxic sites: n=5; normal sites: n=6) and mean residual index of overall body condition, 
calculated as wet weight divided by total length (analysis of variance (ANOVA): P>0.05; hypoxic sites: n=6; normal 
sites: n=10), are provided. Sampling sites served as replicates for all analyses. Error bars indicate the standard errors 
of the means. For (B) female and (C) male spot, box plots present residual GSI values by site and site type. Numerals 
along the x-axis indicate the numbers assigned to sampling sites by the Southeast Area Monitoring and Assessment 
Program and the year during which specimens were collected. The upper and lower parts of each box represent the 
first and third quartiles (the 25th and 75th percentiles), and the horizontal line is the median. The whiskers extend 
above and below the box no more than 1.5 times the interquartile range. Data points inside and outside the boxes 
represent residual GSI values for individual fish. 


the discovery of transitional individuals, highlights the 
need for further research into this phenomenon. Fish spec- 
imens examined in previous work were sampled later in 
the season and were more reproductively developed than 
the fish examined in this study (Thomas and Rahman, 
2012), indicating that hypoxia-induced reproductive 
impairment in fish in the GOM is more widespread and 
develops earlier than previously known. It is important to 
note that, although the results of this study indicate that 
ovarian masculinization in fish can occur earlier in gonadal 
development than previously thought, more research is 


needed to understand how this early impairment trans- 
lates to reproductive success later in the season. 


Ovarian masculinization and reproductive impairment 


For both female and male Atlantic croaker, no statisti- 
cally significant differences were found in the residual 
GSI between hypoxic sites and sites with normal oxy- 
gen concentrations, although a trend of increased GSI 
in specimens from normal sites was observed. This 
finding contradicts those of previous studies in which 
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statistically significant decreases in GSI values were 
found for Atlantic croaker from hypoxic sites (Thomas 
and Rahman, 2012). Although unexpected, the results 
of our study could be explained by the difference in the 
severities of hypoxic conditions between the investiga- 
tions, given that hypoxia in 2016 was less patchy and 
more widespread compared to that in 2007 (the year of 
Thomas and Rahman’s field collections), or by differences 
in the sexual maturity of sampled specimens due to dif- 
ferences in how far into the seasons fish were sampled. 
These discrepancies highlight the importance of contin- 
ued research into hypoxia and its effects, especially those 
on fish at different levels of reproductive maturity. 

Both female and male spot from hypoxic sites had sig- 
nificantly lower (P<0.05) residual GSI values than those 
for females and males from sites with normal oxygen lev- 
els. This finding supports observations from previous stud- 
ies of Atlantic croaker and further highlights the potential 
effects of an environmental stressor, such as hypoxia, on 
fish reproduction at population and ecosystem levels. 

Female and male Atlantic croaker and female spot sam- 
pled from hypoxic sites in 2017 had some of the highest 
GSI values observed for both species and both sexes (Figs. 6 
and 8). The patchy areas of hypoxia from which these spec- 
imens were sampled in 2017 were a fraction of the size of 
the large, more continuous hypoxic zones sampled in 2016 
and 2019 (the hypoxic zone was patchier in 2019 than in 
2016, but the patches of hypoxia were much larger in 2019 
than in 2017) (Fig. 1). The patchiness of hypoxia in 2017 
could have led to hypoxia exposure that was less severe 
than that to which fish were exposed in 2016 and 2019. 
Although we have no diet data to support this interpre- 
tation, it is possible that the smaller areas of hypoxia in 
2017 could have allowed benthic fish to feed on hypoxia- 
stressed infaunal invertebrates, similar to what has been 
observed in the more episodic zones of hypoxia that develop 
in Chesapeake Bay (Pihl et al., 1992). This interpretation 
could also explain why fish from hypoxic sites in 2017 had 
such high GSI values. Balancing the potential for increased 
foraging opportunities associated with hypoxic conditions 
with the physiological costs associated with hypoxia expo- 
sure is one explanation for the context-dependent nature of 
thresholds for hypoxia avoidance in groundfish in the GOM 
(Craig, 2012), highlighting the need for more research to 
place the reproductive costs of hypoxia in the broader eco- 
logical context of the GOM. 


Conclusions 


The results of our study indicate that reproductive impair- 
ment, specifically the presence of intersex individuals 
associated with hypoxic conditions, affects more species 
of fish than previously documented. We also found evi- 
dence of reproductive anomalies, such as the presence of 
ovum remnants in early stages of ovarian development, in 
fish in and around the hypoxic zone in the GOM in 2016. 
These findings, along with the lack of ovarian masculin- 
ization and reproductive anomalies in reference fish and 


the significantly lower GSI values for spot from hypoxic 
sites than from sites with normal oxygen concentrations, 
strongly support the idea that hypoxic conditions can alter 
the reproductive health of fish. 

The findings of this study expand upon previous reports 
of ovarian masculinization and reproductive impairment 
in fish in the GOM. They increase the number of species 
affected and provide novel examples and pathologies of 
this phenomenon. This study also resulted in evidence of 
ovarian masculinization occurring in fish sampled at sites 
bordering hypoxic areas and in fish at earlier developmen- 
tal stages than have been reported previously. However, 
these findings need to be further explored to fully under- 
stand how they fit within the broader context of hypoxia- 
induced ovarian masculinization in fish in the GOM. 


Resumen 


Los resultados de este estudio aportan evidencia adicional 
de masculinizaci6n ovarica en la corvina (Micropogonias 
undulatus) dentro de la zona hipéxica del Golfo de México 
y sus alrededores, asi como las primeras pruebas de mas- 
culinizaciOn ovarica en la croca (Leiostomus xanthurus) y 
el rubio oj6n (Prionotus longispinosus). En algunos ejem- 
plares de corvina y croca observamos este fenédmeno en un 
grado que no se habia reportado antes en ninguna de las 
dos especies. También, encontramos evidencia de génadas 
mas pequenas en crocas de sitios hipéxicos en comparacion 
con las génadas crocas de sitios con concentraciones nor- 
males de oxigeno. Estos hallazgos indican que la masculin- 
izacion ovarica y la deficiencia reproductiva asociados con 
la exposicion a la hipoxia, afectan a una gama de especies 
mas amplia de lo que se conocia hasta ahora y destacan la 
importancia de seguir investigando sobre este tema. 
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Abstract—A bottom-trawl survey with 
a stratified-random sampling design 
has been used to inform stock assess- 
ments for commercially important spe- 
cies in the Gulf of Alaska since 1984. 
A new stratified sampling design was 
evaluated to determine whether its 
use could improve the precision and 
accuracy of abundance estimates. In 
this proposed approach to defining 
strata, historical survey data are used 
to generate what we refer to as infor- 
mation scores (ISes). We compared the 
traditional stratification scheme with 
the new method and 2 other sampling 
designs, using both a design-based esti- 
mator and a model-based estimator 
with each design, to determine if the 
existing approach is optimal. Statistical 
robustness, measured in terms of coef- 
ficient of variation, bias, and root mean 
square error, was compared among 
7 scenarios with different combinations 
of estimators and sampling designs by 
using simulation with a_ spatiotem- 
poral generalized linear mixed model 
conditioned on historical observations 
of catch per unit of effort of 3 species. 
The combination of the design-based 
estimator with the IS-based stratifi- 
cation scheme was the best scenario 
across all performance metrics for all 
species. This scenario consistently had 
the lowest variance and smallest total 
error, and it was generally unbiased. 
In contrast, the pairing of the model- 
based estimator with this sampling 
design was by far the worst-performing 
scenario. The performance of the exist- 
ing approach was average. 
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Surveys for estimating abundance of 
natural resources, and the precision of 
abundance estimates based on survey 
data, are limited in scope by a number 
of logistical constraints, including but 
not limited to funding, staffing, time, 
and access to appropriate and suffi- 
cient sampling tools and platforms. 
The precision of abundance estimates, 
including those of biomass, generally 
increases with sample size, which is 
limited by available resources, assum- 
ing the availability of a species to sam- 
pling remains constant. In practice, the 
scope of a survey effort, as indicated by 
the number of stations where sampling 
occurs, is typically known before the 
sampling design is generated and is 
dictated by logistical constraints. The 
goal of a sampling design, therefore, 
is to optimize the geographical dis- 
tribution of the predetermined num- 
ber of survey stations to achieve the 
highest possible precision for survey 
data products. A potential problem for 
survey continuity arises when circum- 
stances compel a reduction or increase 
in sampling effort, especially after a 
survey time series is well established. 


In such an event, it may be appropri- 
ate to review the existing sampling 
design of the survey for optimality and 
to consider changes to the design while 
assuring continuity of the time series. 
The long-established and relatively 
large-scale multispecies bottom-trawl 
survey that has been conducted in the 
Gulf of Alaska (GOA) since 1984 can 
serve as an illustrative example of how 
these issues can be addressed. 

The NOAA Alaska Fisheries Sci- 
ence Center has conducted a bottom- 
trawl survey in summer in the GOA 
biennially since 1999 (and _ trienni- 
ally between 1984 and 1999), using a 
stratified-random sampling design to 
assess the distribution and abundance 
of groundfish for the purposes of fish- 
eries management (von Szalay and 
Raring, 2018). Traditionally, the survey 
effort consisted of 3 vessels sampling 
approximately 820 stations in 59 differ- 
ent strata. Prior to 1990, the focal sur- 
vey areas in a given year were adapted 
to meet specific scientific or manage- 
ment needs, and as a result the effort 
levels in the different strata could be 
highly variable from one survey year 
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to another. This practice was discontinued in 1990, when 
the sampling design of the survey was standardized to 
adhere to a modified Neyman allocation scheme. A variety 
of factors in recent years, such as funding limitations and 
inability to acquire a sufficient number of suitable char- 
ter vessels, have resulted in the need to reduce the survey 
effort in some years (von Szalay, 2015). Other fisheries 
resource management agencies, both domestic and inter- 
national, similarly have had to or are preparing to contend 
with potential unavoidable reductions in survey effort for 
a variety of reasons, including vessel breakdowns, lack 
of enough survey vessels or staff due to the COVID-19 
pandemic, bad weather, and failure to obtain sampling 
permits for all parts of the survey area (Peel et al., 2013; 
ICES, 2020). 

The reduction of survey effort has amplified the need 
not only to ensure the use of the optimal sampling design 
but also to seek the most robust and precise abundance 
estimator possible that would perform well in years with 
fewer data. This issue can be addressed by comparing the 
statistical robustness of the existing approach with that of 
alternative scenarios based on different combinations of 
estimators and sampling designs, a valuable exercise even 
if survey effort has not been reduced. Until recently, it has 
been difficult to make direct comparisons between differ- 
ent estimators and sampling designs because of a lack of 
an appropriate simulation framework. However, Kotwicki 
and Ono (2019) developed a framework for simulating dis- 
tributions of fish population density with a spatiotemporal 
model conditioned on historical catch and environmental 
survey data for a number of commercially important spe- 
cies in the GOA. We used this operating model to simu- 
late the “true” distribution of abundance and the results of 
simulated sampling with different designs. 

When the existing stratified-random sampling design 
was first devised in the early 1980s, no historical catch or 
environmental data were available to refine the stratifica- 
tion scheme. Now that a long-running time series is well 
established, it is possible to use these data to explore alter- 
native stratification schemes and sampling designs, which 
may potentially reduce the uncertainty in abundance esti- 
mates by yielding more accurate and precise estimates 
of biomass and its associated variance. In recent studies, 
optimization-based methods of stratifying finite popula- 
tions (Hidiroglou and Kozak, 2018; Oyafuso et al., 2021, 
2022) have been proven to be superior to the traditional 
stratification scheme used for the bottom-trawl survey in 
the GOA. In our study, we evaluated the performance of 
3 different traditional sampling designs, which are among 
the most commonly used in the world and have already 
been implemented by the Alaska Fisheries Science Cen- 
ter: simple random, stratified random, and systematic. 

For these evaluations, we used a proposed new technique 
for determining stratum boundaries in the GOA, one in 
which estimates of abundance and its variability from his- 
torical survey data are used to derive an empirical metric 
referred to herein as an information score (IS). This tech- 
nique, which ensures that more stations are allocated to 
strata with higher abundance and variance in abundance 


based on historical survey data, fine-tunes the stratifica- 
tion process by decreasing the variance within strata while 
increasing the variance among strata. This approach con- 
trasts with the traditional stratification scheme whose 
performance may be compromised by including an ad hoc 
management area and depth bounds to define strata. Man- 
agement area boundaries are unlikely to have any bearing 
on fish density, and arbitrary depth boundaries in regular 
100-m intervals are unlikely to coincide with breaks in fish 
density distributions. 

There has recently been a rapid development of meth- 
ods for estimating population density by using spatial 
model-based estimators to calculate biomass or abun- 
dance (e.g., Thorson and Barnett, 2017; Anderson et al.’). 
Design-based estimators, such as the one developed by 
Wakabayashi et al. (1985), are simpler and faster, but spa- 
tial model-based estimators, such as vector autoregres- 
sive spatiotemporal (VAST) models (Thorson, 2019), may 
be more accurate or precise under some scenarios because 
they allow for spatial and temporal correlation and covari- 
ate effects. In addition, model-based designs can be both 
efficient and flexible in allowing uneven sampling due to 
survey logistics and in providing a general framework to 
answer specific design questions (Peel et al., 2013). 

The objective of this study was to compare the statisti- 
cal robustness of 7 scenarios with different combinations 
of sampling designs and estimators in order to identify 
whether an approach improves the accuracy and precision 
of estimates of biomass, which we used as a measure of 
abundance. As a case study, we evaluated the statistical 
robustness of 7 scenarios by simulating the distribution 
and abundance of 3 species in the GOA: Pacific cod (Gadus 
macrocephalus), Pacific ocean perch (Sebastes alutus), and 
arrowtooth flounder (Atheresthes stomias). These species 
were chosen on the basis of their ecological or commercial 
importance as well as the diversity of their spatial distri- 
bution patterns. In 6 scenarios, a random sampling design 
was paired with a design-based or model-based estimator, 
and in 1 scenario a systematic sampling design was paired 
with a model-based estimator. 


Materials and methods 
Survey characteristics 


The GOA forms the northeastern border of the Pacific 
Ocean and has complex bathymetric features ranging 
from jagged, mountainous pinnacles to flat, muddy areas 
(von Szalay and Raring, 2018). These features provide a 
variety of habitats, creating a complex ecosystem. The 
standard survey area of the biennial bottom-trawl survey 
in the GOA is approximately 320,000 km” and includes 


1 Anderson, S. C., E. J. Ward, P. A. English, and L. A. K. Barnett. 
Unpubl. manuscript. sdmTMB: an R package for fast, flexi- 
ble, and user-friendly generalized linear mixed effects mod- 
els with spatial and spatiotemporal random fields. bioRxiv 
2022.03.24.485545. [Available from website. ] 
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Figure 1 


Map showing the 59 strata of the bottom-trawl survey conducted biennially in the Gulf of Alaska 
by the NOAA Alaska Fisheries Science Center. The survey area includes the portion of the con- 
tinental shelf from the Islands of Four Mountains eastward approximately 2800 km to Dixon 
Entrance and from nearshore waters to a depth of 1000 m. The stratum boundaries have been 


fixed since the inception of the survey in 1984. 


the portion of the continental shelf from the Islands of 
Four Mountains, in the Aleutian Islands, eastward to 
Dixon Entrance, near the border of Canada and the United 
States, and from nearshore waters to a depth of 1000 m 
(Fig. 1). The continental shelf, which makes up 82% of the 
survey area, is bordered by the continental slope, a region 
that is approximately 20 km in width, but only the portion 
of the slope at depths between 200 and 1000 m, is sampled 
in the standard surveys. The total survey area has been 
effectively decreased to 308,415 km? in recent surveys as 
reductions in survey effort have resulted in the deepest 
strata (700-1000 m) not being sampled. 

In the bottom-trawl surveys conducted in the GOA by 
the Resource Assessment and Conservation Engineering 
Division of the Alaska Fisheries Science Center, approx- 
imately 820 locations within 59 strata traditionally have 
been sampled with a stratified-random sampling design. 
Abundance estimates are generated by using the estima- 
tor developed by Wakabayashi et al. (1985). Henceforth, 
this combination of sampling design and estimator will be 
referred to as the TRS. 

The same standard trawl gear has been used since 
the beginning of the survey, a Poly Nor’Kastern 4-seam 
bottom trawl with 24.2-m roller gear constructed with 
36-cm rubber bobbins separated by 10-cm rubber disks 
(Stauffer, 2004). Surveys start in the western end of 


the survey area and proceed eastward. Tow duration is 
approximately 15 min at 1.54 m/s (8 kt). The catch per 
unit of effort (CPUE) is estimated by using the area-swept 
method (e.g., Alverson and Pereyra, 1969), which defines 
the effort as the product of the distance fished and the 
average distance between wing tips (for details, see von 
Szalay and Raring, 2018). 

A grid with a resolution of 5 km superimposed on the 
survey area is used for sampling design. Each 5-by-5-km 
cell in the grid is a potential survey station, but some 
potential stations are excluded from the pool of stations 
eligible for survey selection because they are untrawlable 
as a result of rough bottom or other factors. Further- 
more, because they are truncated by the coastline and the 
deepest edge of the sampling domain, grid cells along the 
boundaries of the survey area have an area smaller than 
a standard 25-km’ grid cell. Only grid cells large enough 
to accommodate a complete tow path (those with an area 
>5 km”) are eligible for sampling. 


Simulated distributions of catch per unit of effort 


Simulated CPUE distributions of Pacific cod, Pacific 
ocean perch, and arrowtooth flounder were created by fit- 
ting a spatiotemporal delta model to historical GOA sur- 
vey data from 1996 through 2015, with the package 
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R-INLA (vers. 18.07.12; available from website, accessed 
October 2018; Rue et al., 2009) in R (vers. 3.6.1; R Core 
Team, 2019); this package accounts for both environmen- 
tal covariates and spatiotemporal dependency in catches 
(Kotwicki and Ono, 2019) (Fig. 2). The delta model has 
2 components: one that models the species occurrence 
and another that models positive CPUE. 

Species occurrence at locations s during year f¢, 1,(s), 
was modeled by using a binomial generalized linear mixed 
model with the logit link function (see Lindgren et al., 2011): 


logit (m, (s)} =X, (s)b + Wy (s), where (1) 


O, (s) ~ N (p10, 6a 


and X,(s) = the matrix of covariates at locations s during 
year J; 
6 = the vector of regression coefficients; 
w= the spatiotemporal variation that follows a 
first-order autoregressive process; 
N = normal distribution; 
0, = the degree of autoregression in encounter prob- 
ability between successive years; and 
x, = the spatial covariance modeled as a Matern 
function with smoothness of 1. 


58°N aumes ———___—- Bs as Sees 
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The covariates used in this study were log(depth), 
(log(depth))*, bottom temperature, bottom temperature 
squared, surface temperature, surface temperature squared, 
and a fixed year effect (Kotwicki et al., 2005, 2015). 

The non-zero species density at a set of locations s 
during year ft, p,(s), was modeled by using a lognormal 
distribution: 


log (44 (s))= Zi, (s)a + 6,(s), where (2) 


o; (s) ~ N (p23, (s), 25] 


and Z,(s) =a matrix of covariates at locations s during 
year f; 
a = the vector of regression coefficients; 

6, = the spatial field for year t assumed to follow an 
autoregressive 1 process; 

0, = autocorrelation of the autoregressive 1 pro- 
cess in which the current value is based on the 
immediately preceding value; and 

x, = the spatial covariance modeled as a Matern 
function with smoothness of 2. 


Annual distributions of simulated CPUE within the 
survey area were predicted over a nominal grid with a 


CPUE (kg/km?) 
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Figure 2 


Simulated distribution of catch per unit of effort (CPUE) of arrowtooth flounder (Atheresthes sto- 
mias) in 2015 at the entrance to Cook Inlet in southcentral Alaska. The simulations were created 
by fitting a spatiotemporal delta model that accounts for both environmental covariates and spa- 
tiotemporal dependency in catches to historical data from the bottom-trawl survey in the Gulf 
of Alaska. Each pixel on the map is 2 km” and represents the CPUE realized when a station is 
sampled at that location in the simulated surveys conducted in this study. 
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2-km resolution on the basis of a single but random sam- 
ple from the posterior predictive distribution of the mod- 
els while accounting for the sampling process. In other 
words, simulated CPUE was calculated as the product of 
sampled fish occurrence and sampled density. Sampled 
fish occurrence was estimated by using a Bernoulli trial 
with the probability determined by the parameter sam- 
ples taken from the posterior predictive distribution Gin 
practice, a random Markov chain Monte Carlo [MCMC] 
run was chosen and parameter values were taken from 
it). Sampled density was estimated by using a sample 
from a Gaussian distribution on the link scale with the 
mean and variance derived from parameter values taken 
from a randomly chosen MCMC run (the sample value 
was exponentiated by using the natural exponential func- 
tion to convert it to real space). 

With this approach, predictions accounted for all 
sources of uncertainty included in the model and created 
a noisier CPUE distribution that was more reflective of 
“true” patterns than the mean MCMC prediction. Envi- 
ronmental covariate values were determined by kriging 
with the semivariogram model that best fit the historical 
survey data (Kotwicki and Ono, 2019), as implemented by 
the function autofitVariogram in the R package automap 
(vers. 1.0-16; Hiemstra et al., 2009). All geographic coordi- 
nates were converted into an Albers projection in order to 
preserve distances prior to analysis. 


Simulation of surveys 


The CPUE distributions defined over the grid with the 
finer resolution of 2 km were superimposed on the grid 
with a 5-km resolution for the surveys in the GOA before 
the simulated surveys were generated. Each station sam- 
pled during surveys in the GOA (i.e., each 25-km? grid 
cell), therefore, contained an average of approximately 
4—7 of the 2-km’ grid cells with density data. A variety of 
scenarios with different combinations of estimators and 
sampling designs were then implemented to generate sim- 
ulated surveys over the modeled CPUE distributions for 
each of the 10 survey years between 1996 and 2015 and 
for each of 3 commercially important species: Pacific cod, 
Pacific ocean perch, and arrowtooth flounder. 

In addition to the traditional stratified-random sam- 
pling design and the proposed stratification sampling 
design based on an IS, we considered 2 additional sam- 
pling designs, simple random sampling and systematic 
sampling, in conjunction with both the model-based 
estimator and the design-based estimator, for a total of 
8 potential scenarios. The scenario in which the design- 
based estimator is paired with the systematic sampling 
design was excluded from this analysis because of its 
unsuitability for the survey area in the GOA. Because the 
survey area is characterized by relatively wide strata on 
the continental shelf but extremely narrow strata along 
the continental slope, stations are too far apart to ade- 
quately cover the slope area. The remaining 7 scenarios 
consisted of the use of the design-based estimator with 
3 different sampling designs—simple random sampling 


(SRS), stratified-random sampling based on ISes (WIS) 
(described later in the “Stratified-random sampling based 
on information scores” section), and TRS—and the use of 
the model-based estimator in the R package VAST (vers. 
3.3.0; Thorson, 2019) with 4 different sampling designs— 
systematic sampling (VSY), simple random sampling 
(VRS), stratified-random sampling based on ISes (VIS), 
and traditional stratified-random sampling (VTRS). 

We simulated 100 replicate surveys, each consisting of 
820 stations, for each of the 10 years in which surveys had 
been conducted, resulting in 82,000 simulated stations per 
survey year for each scenario. All analyses were conducted 
in R, vers. 3.6.1 (R Core Team, 2019). 


Design-based approaches 


Traditional stratified-random sampling In the stratified- 
random sampling design used since the inception of the 
bottom-trawl survey of the GOA in 1984, the survey area 
is divided into 59 strata defined by water depth, bottom 
terrain (e.g., shelf, gully, and slope), and statistical man- 
agement areas (von Szalay and Raring, 2018). Following 
this traditional design, we allocated stations among the 
strata for species k by using a modified Neyman strategy 
for optimal allocation (Cochran, 1977): 
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where n;;, = the sample size for species k in stratum h; 
n = the total sample size; 
N,, = the population size for stratum h; 
Sy, =the standard deviation of species k in 
stratum h; 
a,, = the area of stratum h; and 
Cc), = the cost to conduct a trawl tow in stratum h. 


hy = 


Catch rates, stratum variances, and stratum areas from 
the surveys conducted in 1990-2015 were used to allocate 
sampling effort among strata for each previous survey 
year and for each of the 50 principal GOA groundfish spe- 
cies that the historical survey was designed to sample. The 
estimated time to perform a trawl tow in a given stratum 
was used as a cost variable because tows in deeper strata 
have a greater probability of unacceptable gear perfor- 
mance and take longer to complete. A mean sample size 
was calculated across years for each principal species, and 
the sample size for each species was then weighted by each 
species’ commercial value, which was defined as the prod- 
uct of the mean biomass and its ex-vessel value in 2015. 
The calculated sample sizes, representing the numbers of 
trawl tows allocated to the various strata, were rounded 
to whole numbers while ensuring that each stratum was 
allocated at least 2 samples. 

The allocated stations within each stratum were 
randomly selected without replacement from the 
2-km-resolution grid described earlier. Sample mean (*},) 
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and standard deviation (S},,) for species k were calculated 
by stratum and were used to generate the survey mean 


(Xr) and variance (S35) for a stratified-random sampling 
design according to the following equations: 
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Stratified-random sampling based on information scores As 
an alternative method of stratification to the traditional 
scheme, we propose the IS approach to stratified-random 
sampling. This novel approach uses spatial and temporal 
information of abundance and temporal information of 
variation in abundance to derive an IS for each sampling 
unit. A high IS should identify areas where abundance 
tends to be high, areas where the abundance is highly 
variable over time, or a combination of these characteris- 
tics. Sampling in such areas is likely to provide most infor- 
mation about population abundance, composition, and 
trends in abundance, all of which are critical for use in 
stock assessment models and for acquisition of other types 
of information about sampled populations and temporal 
trends in population parameters. Information scores could 
be used to stratify the survey area and provide a basis for 
prioritizing areas with high fish density and high variance 
in fish density to ensure that more stations are allocated 
to strata with higher expected abundance and variability 
in abundance. Given this logic, we propose an IS that is 
calculated for each of the 15,718 eligible 5-by-5-km grid 
cells in the survey area as follows: 
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where CPUE;, = the mean predicted CPUE of species k 
in grid cell 1 over the 10 survey years 
between 1996 and 2015. 


The CPUE\, of grid cell i (5-by-5 km) for a given year was 
calculated as the mean prediction of the simulated densi- 
ties of all 2-by-2-km grid cells within it. By squaring the 
mean predicted CPUE, we ensure that the 2 terms have 
the same units and equal weight. 

To determine how to stratify on the basis of discretiz- 
ing the ISes into bins of similar values, we first examined 
how the scores were distributed. We found the distribu- 
tion of the 15,718 ISes to be highly positively skewed 
for all species (Pearson’s coefficient of skewness ranged 
between 8.7 for arrowtooth flounder and 43.5 for Pacific 
ocean perch; Suppl. Fig. 1 [online only)). 

These values were sorted from low to high, and evenly 
spaced breakpoints were inserted to generate 20 strata, 
each with 786 stations. The choice to use 20 strata resulted 
from observations about the scores. These observations 


were that approximately 5% of the ISes were equal to zero, 
and a similar proportion of the scores on the opposite end 
of the distribution were of the same order of magnitude, 
resulting in 2 strata that contained the extreme values 
on both ends of the distribution with other strata having 
intermediate values. In a sensitivity analysis, we found no 
indication of substantial differences in results when either 
19 or 21 strata were used. Furthermore, 20 strata is a rea- 
sonable number for the survey in the GOA and is also the 
number of strata being considered for a revised sampling 
design for this survey. 

The number of survey stations allocated to the different 
strata was weighted by the mean IS of each stratum. First, 
we used the Neyman algorithm to determine the sample 
allocation (Lavrakas, 2008): 
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Then, to weight the allocation among strata by using the 
mean IS of each stratum, the weights were calculated as 


follows: 
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where w), = the weight of stratum h, and 


IS; = the mean information score for stratum h. 


In this equation, substituting J 1Sn for s, in Equation 
7 and algebraically simplifying the new equation (N,, is 
a constant and cancels out) results in n,, being equal to 
nwy,, which is the sample size of stratum h weighted by its 
information score. 

Some stations originally assigned to the strata with the 
highest ISes were reallocated to strata receiving fewer 
than 2 stations to make variance calculations possible 
for all strata. This procedure for station allocation among 
strata was repeated for each species. 

The allocated stations within each stratum were ran- 
domly selected without replacement from the grid with a 
5-km resolution for the survey in the GOA. A random grid 
cell from the 2-km-resolution grid was then sampled to 
represent the CPUE of the entire station. Sample mean 
and standard deviation were calculated by stratum and 
were used to generate the survey mean and variance for 
a stratified-random sampling design similar to the proce- 
dure used for the traditional stratified-random approach 
described in the previous section. 


Model-based approaches 


The model-based estimator was a spatiotemporal delta 
model consisting of a binomial model for the probability 
of encounter and a user-selected model for positive catch 
rates, as implemented in the VAST package (Thorson, 
2019). In a report of their recent work, Thorson et al. (2021) 
noted that the choice of distributional assumptions can 
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have a substantial effect on the scale of the biomass index 
estimated from such spatiotemporal models; therefore, we 
tested a range of assumptions to determine which provided 
the best fit and accuracy. We considered 2 different obser- 
vation error distributions for positive catch rates, lognor- 
mal and gamma, and whether or not to correct values for 
retransformation bias (Thorson and Kristensen, 2016). 

The hurdle models included year as a fixed effect, a 
spatial random field for each model component (encoun- 
ter and positive catch rate), and independent spatiotem- 
poral random fields for the positive catch rate component 
(Suppl. Table) (online only). The spatiotemporal term for 
encounter probability was not estimated because of the 
lack of consistent convergence among replicate simulated 
surveys. A model resolution of 500 knots was used, and 
bilinear interpolation was implemented to extrapolate 
from knot locations in order to generate predictions at 
each location in the same 2-km-resolution grid used in the 
operating model. Anisotropy was estimated. The models 
for each scenario differed only in terms of the observation 
error model used for positive catch rates (gamma versus 
lognormal distribution) and by whether retransformation 
bias was addressed (i.e., with or without epsilon bias cor- 
rection) (Thorson and Kristensen, 2016). All models were 
run with the VAST package (Thorson et al., 2015) and the 
tools available in FishStatsUtils (vers. 2.6.0; available 
from website, accessed October 2019). 

A customized 2-km-resolution extrapolation grid was 
created as a subset of the standard GOA survey grid 
(depths <1000 m), which excluded all grid cells known to 
be untrawlable. This custom grid was used to ensure that 
the VAST estimation model only integrated densities from 
trawlable areas or areas with an unknown trawlability 
status in order to make the estimates of biomass and its 
variance from the model-based approaches comparable 
to those from the design-based approaches. The results 
of the 4 different combinations of settings, between the 
2 observation models for positive catch rates and whether 
or not the retransformation bias was addressed, for each 
scenario and species were compared with respect to the 
performance metrics described later in the “Performance 
metrics” section. The combination of settings that resulted 
in the lowest relative root-mean-square error (rRMSE) for 
a given scenario and species was selected for final compar- 
ison with the other scenarios (Table 1). 


Simple random sampling Simulated CPUE values from 
the operating model were sampled by using simple ran- 
dom sampling. These values were then used as input to 
the VAST estimation model, along with the associated 
data on location and year for CPUE observations. 


Systematic sampling A rectangular grid encompassing 
the entire survey area in the GOA was created by identi- 
fying the extreme latitudes and longitudes spanning the 
survey area. An iterative procedure was used to deter- 
mine that the length of the square grid cells that would 
result in approximately 820 uniformly spaced stations 
falling within the survey area was approximately 20 km. 


To generate multiple iterations of systematic surveys, the 
origin of the survey grid was altered between iterations by 
randomly sampling a set of coordinates from within the 
erid cell in the extreme southeast corner of the survey grid. 
This 20-km? survey grid was, in turn, superimposed on the 
nominal 2-km-resolution grid of simulated fish densities. 
The 2-km? density grid cells within each of the 820 uni- 
formly spaced grid cells were identified, and the density 
associated with the 2-km’ grid cell closest to the center of 
the survey grid cell was used to represent the density of 
the entire survey grid cell. 


Stratified-random sampling based on information scores 
Simulated CPUE estimates from the operating model 
were sampled by using an IS-based stratified-random 
sampling design. These values were then used as input 
to the VAST estimation model, along with the associated 
data on location and year for CPUE observations. 


Traditional stratified-random sampling Based on the same 
principles described for the traditional stratified-random 
sampling in the “Design-based approaches” section, the sta- 
tion allocation by stratum scheme replicated the sampling 
design used for the bottom-traw] survey in the GOA in 2007 
that consisted of 825 stations allocated to the 59 survey 
strata. However, rather than sampling directly from the 
2-km-resolution CPUE grid, the stations were sampled from 
the grid with a 5-km resolution used for the survey in the 
GOA, and the CPUE associated with each station was based 
on the mean of the 2-km? CPUE grid cells within it. 


Performance metrics 


Three performance metrics were used to evaluate the rel- 
ative performance of the 7 scenarios that were used to 
simulate abundance: coefficient of variation (CV), relative 
bias, and rRMSE. Each metric was applied to estimates 
of biomass and its associated variance and was calculated 
separately for each year and species. The true biomass was 
calculated as the product of the trawlable area of the sur- 
vey area in the GOA and the arithmetic mean of simulated 
CPUE (with observation error) from all 65,863 trawlable 
erid cells. The true variance was calculated for each year as 
the variance of the replicate simulated biomass estimates: 


Diy =Bay 
N 


On.= (9) 


9 


where 0 a = a “true” variance of the abundance index based 
on 100 simulated surveys; 
N = the total number of simulated surveys (i.e., 100); 
Bg = the estimated biomass realized in simulation 
survey s; and 
By = the “true” biomass estimated from simulated 
density grids. 


Bias of both biomass and its variance was defined as the 
mean of the deviations between the respective value of 
each simulation and the true value. The relative bias (RB) 
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Table 1 


Relative root-mean-square errors (rRMSEs) for 7 scenarios with different combinations of estimators and 
sampling designs used to estimate biomass and its variance for arrowtooth flounder (Atheresthes stomias) 
(ATF), Pacific cod (Gadus macrocephalus) (COD), and Pacific ocean perch (Sebastes alutus) (POP) in the 
Gulf of Alaska. The observation error model for positive catch rates and the approach to bias correction 
are indicated for the scenarios in which the model-based estimator was used. The model-based estimator 
was the vector autoregressive spatiotemporal (VAST) model in the R package VAST. The design-based esti- 
mator was developed by Wakabayashi et al. (1985). One of the sampling designs was a stratified-random 
design based on information scores. The scenarios with the best overall performance for each species have 
the lowest values for biomass and its variance, which are indicated with asterisks. Analysis of scenarios 
was based on data from the bottom-trawl surveys conducted in the Gulf of Alaska during 1996-2015 by 
the NOAA Alaska Fisheries Science Center. The abbreviations in parentheses represent the 7 scenarios. 


rRMSE 
Observation Bias 


error model correction 


Estimator and sampling 


Species design (scenario) Biomass Variance 

Model-based approaches 

ATF VAST systematic (VSY) 
VAST simple random (VRS) 
VAST information score (VIS) 


VAST stratified random (VTRS) 


VAST systematic (VSY) 

VAST simple random (VRS) 
VAST information score (VIS) 
VAST stratified random (VTRS) 


VAST systematic (VSY) 

VAST simple random (VRS) 
VAST information score (VIS) 
VAST stratified random (VTRS) 


Design-based approaches 

ATF Simple random (SRS) 0.19 
Information score (WIS) 0.05* 
Stratified random (TRS) 0.28 


Simple random (SRS) 0.15 
Information score (WIS) 0.09* 
Stratified random (TRS) 0.19 


Simple random (SRS) 0.39 
Information score (WIS) 0.10* 
Stratified random (TRS) 0.40 


delta-gamma Off >1.00 
delta-lognormal Off 0.30 

delta-gamma On 0.25 

delta-gamma On 0.18* 


delta-lognormal On 0.69 

delta-lognormal On 0.18 
delta-gamma On 0.45" 
delta-gamma On 0.15* 


delta-gamma On 0.98 
delta-lognormal On 0.31 

delta-gamma Off 0.18* 

delta-gamma On 0.28 


of these estimates was calculated within each year and 
species as follows: 


a ~ (8: a Ore 
RB = 2—_____, (10) 
Ore 
where 9; = the value (biomass or variance) of the ith sim- 
ulation replicate, and 
Ore = the true value (biomass or variance), over 100 
simulation replicates. 


The CV was defined as follows: 


vie esl where (11) 


Orrue 


2 
n 


1 
Var(®) = pat 7 Boon | 


and Var(8) = the mean of the square deviations between 
the value (biomass or variance) of each simu- 
lation and the mean biomass or variance. 


The mean square error was defined as the sum of the bias 
squared and variance, and the rRMSE was calculated as 


follows: 
\ Bias? + Var 


rRMSE Se, (12) 


True 


where Bias = the absolute bias as defined by the numera- 
tor in Equation 10. 
Results 


Of the 7 scenarios with different combinations of estima- 
tors and sampling designs that were considered in this 
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study, the best scenario for estimating biomass and its 
variance across all performance metrics and all species 
was the design-based estimator coupled with the stratified- 
random sampling design in which strata were defined by 
ISes. This IS scenario consistently had the lowest CV and 
rRMSE and was generally unbiased (Table 1, Figs. 3-8). In 
contrast, the scenario that paired the model-based estima- 
tor of biomass with this sampling design (VIS scenario) 
was by far the worst-performing scenario by all 3 perfor- 
mance metrics for all species. It had the highest CV and 
rRMSE and was the most biased for all species. However, 
for variance, the VIS scenario performed considerably bet- 
ter, with no bias for 2 of the species and with an rRMSE 
similar to or smaller than that of the other model-based 
variance scenarios. 

The performance of the combination of the design-based 
estimator with the traditional stratified-random sampling 
design (TRS scenario) was generally average for estimates 
of biomass and its variance among the 7 scenarios. The 
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performance of this scenario was most similar to that of 
the scenario that pairs the design-based estimator with a 
simple random sampling design (Figs. 3-8). 

The performance of scenarios differed among species. 
Of the 7 scenarios for estimating biomass, 5 scenarios per- 
formed relatively well for arrowtooth flounder (VRS, VTRS, 
SRS, WIS, and TRS). A sixth scenario (VSY) performed well 
for Pacific cod in terms of rRMSE (rRMSE <0.10). However, 
only the WIS had an rRMSE below 0.10 for Pacific ocean 
perch (Table 1). All of the scenarios with the designed-based 
estimator of biomass were consistently unbiased (mean 
relative bias <0.005 for all species), whereas the scenarios 
with the model-based estimator were generally negatively 
biased, except the VRS for arrowtooth flounder and the VIS 
across all species (Figs. 3-5). 

Of the 7 scenarios for estimating variance in biomass, 
only the WIS consistently performed well as measured by 
rRMSE (rRMSE: 0.06—0.10). Although most of the other 
scenarios were relatively unbiased for Pacific cod and 
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Figure 3 


Relative performance of the 7 scenarios with different combinations of estimators and sampling 
designs in estimating biomass for arrowtooth flounder (Atheresthes stomias) based on 100 simu- 
lated bottom-trawl surveys in the Gulf of Alaska. Each box plot represents the distribution of a 
performance metric, (A) coefficient of variation, (B) relative bias, or (C) relative root-mean-square 
error (RMSE), for a scenario over the 10 survey years between 1996 and 2015. The model-based 
estimator, the vector autoregressive spatiotemporal (VAST) model in the R package VAST, was 
paired with 4 sampling designs: systematic (VSY scenario), simple random (VRS scenario), strati- 
fied random based on information scores (VIS scenario), and traditional stratified random (VTRS 
scenario). The designed-based estimator was paired with 3 sampling designs: simple random (SRS 
scenario), stratified random based on information scores (WIS scenario), and traditional stratified 
random based on Neyman allocation (TRS scenario). In each box plot, the thick black line is the 
median. The lower and upper hinges of the box plots correspond to the 25th and 75th percentiles. 
The whiskers extend from the hinge to the farthest value no more than 1.5 times the interquartile 
range (distance between the 25th and 75th percentiles). Data points beyond the end of the whis- 


kers, outliers, are plotted individually. 
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Pacific ocean perch (with VSY being the only exception), 
they produced relatively large variances, which resulted 
in relatively large rRMSEs (Figs. 6-8). 

The VIRS scenario performed better in estimating bio- 
mass than 2 of the design-based methods (TRS and SRS), 
as well as all of the other model-based scenarios in terms 
of CV and rRMSE, but this scenario was somewhat nega- 
tively biased (mean relative bias ranged between —0.02 
and —0.05). The only scenario with the model-based esti- 
mator of biomass that had a comparable rRMSE was the 
one with the simple random sampling design (VRS sce- 
nario), but only for arrowtooth flounder and Pacific ocean 
perch (Figs. 3-5). The performance of the VTRS in esti- 
mating variance in biomass, on the other hand, was gener- 
ally similar to that of the other methods (except for the 
poorly performing VSY and the best performing WIS) 
across all performance metrics (Figs. 6-8). 

The VSY scenario for estimating the variance of bio- 
mass had extreme values of all performance metrics 
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for Pacific cod and Pacific ocean perch (CV: 0.64—0.78; 
mean relative bias: from —0.28 to —0.60; rRMSE: 0.7—1.0; 
Figs. 7-8) but even more extreme values for arrowtooth 
flounder (CV and rRMSE >1.0). The VSY scenario for esti- 
mating biomass variance for arrowtooth flounder, there- 
fore, was excluded from the boxplots in Figure 6 because 
it was an extreme outlier, which if included, would have 
made comparisons of the other scenarios difficult. The poor 
performance of the VSY scenario for estimating variance 
was due to the highly skewed variance distributions for all 
species, with ranges spanning as much as 5 orders of mag- 
nitude for arrowtooth flounder (Suppl. Fig. 2) (online only). 

Apart from the differences in the relative performance 
of the various scenarios for estimating biomass and vari- 
ance, differences in the performance of each scenario 
among species are notable. All of the scenarios performed 
considerably worse in estimating biomass for Pacific ocean 
perch than for the other 2 species in terms of rRMSE. The 
difference in performance across species was small for the 
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Figure 4 


Relative performance of the 7 scenarios with different combinations of estimators and sampling 
designs in estimating biomass for Pacific cod (Gadus macrocephalus) based on 100 simulated 
bottom-trawl surveys in the Gulf of Alaska. Each box plot represents the distribution of a per- 
formance metric, (A) coefficient of variation, (B) relative bias, or (C) relative root-mean-square 
error (RMSE), for a scenario over the 10 survey years between 1996 and 2015. The model-based 
estimator, the vector autoregressive spatiotemporal (VAST) model in the R package VAST, was 
paired with 4 sampling designs: systematic (VSY scenario), simple random (VRS scenario), strati- 
fied random based on information scores (VIS scenario), and traditional stratified random (VTRS 
scenario). The designed-based estimator was paired with 3 sampling designs: simple random (SRS 
scenario), stratified random based on information scores (WIS scenario), and traditional stratified 
random based on Neyman allocation (TRS scenario). In each box plot, the thick black line is the 
median. The lower and upper hinges of the box plots correspond to the 25th and 75th percentiles. 
The whiskers extend from the hinge to the farthest value no more than 1.5 times the interquartile 
range (distance between the 25th and 75th percentiles). Data points beyond the end of the whis- 


kers, outliers, are plotted individually. 
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Relative performance of the 7 scenarios with different combinations of estimators and sampling 
designs in estimating biomass for Pacific ocean perch (Sebastes alutus) based on 100 simulated 
bottom-trawl surveys in the Gulf of Alaska. Each box plot represents the distribution of a per- 
formance metric, (A) coefficient of variation, (B) relative bias, or (C) relative root-mean-square 
error (RMSE), for a scenario over the 10 survey years between 1996 and 2015. The model-based 
estimator, the vector autoregressive spatiotemporal (VAST) model in the R package VAST, was 
paired with 4 sampling designs: systematic (VSY scenario), simple random (VRS scenario), strati- 
fied random based on information scores (VIS scenario), and traditional stratified random (VTRS 
scenario). The designed-based estimator was paired with 3 sampling designs: simple random (SRS 
scenario), stratified random based on information scores (WIS scenario), and traditional stratified 
random based on Neyman allocation (TRS scenario). In each box plot, the thick black line is the 
median. The lower and upper hinges of the box plots correspond to the 25th and 75th percentiles. 
The whiskers extend from the hinge to the farthest value no more than 1.5 times the interquartile 
range (distance between the 25th and 75th percentiles). Data points beyond the end of the whis- 


kers, outliers, are plotted individually. 


top-performing WIS scenario (an rRMSE of 0.03 for Pacific 
ocean perch versus values of 0.02 and 0.025 for arrowtooth 
flounder and Pacific cod, respectively). However, the differ- 
ence was much higher for all of the other scenarios for esti- 
mating biomass, with absolute differences in rRMSE 
between Pacific ocean perch and the other 2 species rang- 
ing from 0.01 to 0.16 for arrowtooth flounder and from 
0.04 to 0.27 for Pacific cod (Table1). In contrast, all but one 
of the scenarios for estimating variance of biomass per- 
formed best for Pacific cod, with absolute differences in the 
rRMSE between Pacific cod and the other 2 species rang- 
ing from 0.03 to >0.31 for arrowtooth flounder and from 
0.01 to 0.29 for Pacific ocean perch. The single exception 
was the WIS scenario, which had a smaller rRMSE for 
arrowtooth flounder (0.05 for arrowtooth flounder versus 
a value of 0.09 for Pacific cod) but a larger rRMSE for 
Pacific ocean perch (0.10). The differences in relative per- 
formance of the scenarios between Pacific ocean perch and 
arrowtooth flounder were mixed. 


Discussion 


The objective of this study was to assess the performance of 
7 scenarios with different combinations of sampling designs 
and estimators and to compare their statistical robustness. 
Although we used 3 different performance metrics to eval- 
uate statistical robustness, for conciseness, we highlight 
herein the best-performing scenarios on the basis of only 
rRMSE (those with the lowest rRMSEs are considered the 
best performing), which includes both the bias and variance 
components teased apart by the other performance metrics. 

The results of the simulation analyses indicate that, 
compared to the scenario that paired the design-based 
estimator with the existing sampling design (TRS), 4 sce- 
narios performed better for arrowtooth flounder, 3 sce- 
narios performed the same or better for Pacific cod, and 
3 scenarios performed better for Pacific ocean perch. Of 
these scenarios, only the design-based WIS and the 
model-based VTRS scenarios performed better than the 
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Figure 6 


Relative performance of the 7 scenarios with different combinations of estimators and sampling 
designs in estimating the variance of biomass for arrowtooth flounder (Atheresthes stomias) based 
on 100 simulated bottom-trawl surveys in the Gulf of Alaska. Each box plot represents the distri- 
bution of a performance metric, (A) coefficient of variation, (B) relative bias, or (C) relative root- 
mean-square error (RMSE), for a scenario over the 10 survey years between 1996 and 2015. The 
box plots for the VSY scenario for this species have been omitted for readability of the vertical axis 
values of the other methods. The model-based estimator, the vector autoregressive spatiotemporal 
(VAST) model in the R package VAST, was paired with 4 sampling designs: systematic (VSY sce- 
nario), simple random (VRS scenario), stratified random based on information scores (VIS scenario), 
and traditional stratified random (VTRS scenario). The designed-based estimator was paired with 
3 sampling designs: simple random (SRS scenario), stratified random based on information scores 
(WIS scenario), and traditional stratified random based on Neyman allocation (TRS scenario). In 
each box plot, the thick black line is the median. The lower and upper hinges of the box plots corre- 
spond to the 25th and 75th percentiles. The whiskers extend from the hinge to the farthest value no 
more than 1.5 times the interquartile range (distance between the 25th and 75th percentiles). Data 
points beyond the end of the whiskers, outliers, are plotted individually. 


TRS scenario for all 3 species, with only the WIS scenario 
consistently performing substantially better than the 
TRS scenario for estimates of both biomass and its vari- 
ance for all species (with rRMSEs >50% lower than those 
for the TRS). The reduction in rRMSE from the WIS sce- 
nario to the TRS scenario ranged among species between 
58% and 81% for biomass and between 53% and 82% for 
variance in biomass, with the greatest reduction observed 
for Pacific ocean perch. It is likely that the traditional 
sampling design is not as suitable for assessing abun- 
dance of Pacific ocean perch because of the patchy distri- 
bution of this species in the GOA (Lunsford et al., 2001; 
Clausen and Fujioka, 2007), a pattern that requires 
larger sample sizes in areas of high density than in other 
areas. In summary, our results indicate that the combi- 
nation of stratified-random sampling with model-based 
estimation methods is likely to produce the best results. 


These results have at least 2 important implications for 
the redesign of and improvements to the method of abun- 
dance estimation for the bottom-trawl survey conducted in 
the GOA. First, in the short term, estimates of abundance 
can be improved immediately by using model-based estima- 
tors in combination with traditional stratified-random sam- 
pling. Second, in the long term, basing stratification on [Ses 
could be an additional change that can improve estimates 
for species with spatially restricted or patchy distributions, 
when used in combination with design-based estimators. 


Recommendation to develop model-based estimators 


Our results indicate clear improvements in estimation for 
all species when the model-based estimator was used with 
a stratified-random sampling design. This finding is 
important because it indicates that some of the losses in 
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Figure 7 


Relative performance of the 7 scenarios with different combinations of estimators and sampling 
designs in estimating the variance of biomass for Pacific cod (Gadus macrocephaluls) based on 
100 simulated bottom-trawl surveys in the Gulf of Alaska. Each box plot represents the distri- 
bution of a performance metric, (A) coefficient of variation, (B) relative bias, or (C) relative root- 
mean-square error (RMSE), for a scenario over the 10 survey years between 1996 and 2015. The 
model-based estimator, the vector autoregressive spatiotemporal (VAST) model in the R package 
VAST, was paired with 4 sampling designs: systematic (VSY scenario), simple random (VRS sce- 
nario), stratified random based on information scores (VIS scenario), and traditional stratified 
random (VTRS scenario). The designed-based estimator was paired with 3 sampling designs: sim- 
ple random (SRS scenario), stratified random based on information scores (WIS scenario), and tra- 
ditional stratified random based on Neyman allocation (TRS scenario). In each box plot, the thick 
black line is the median. The lower and upper hinges of the box plots correspond to the 25th and 
75th percentiles. The whiskers extend from the hinge to the farthest value no more than 1.5 times 
the interquartile range (distance between the 25th and 75th percentiles). Data points beyond the 


end of the whiskers, outliers, are plotted individually. 


precision expected from design-based estimators in cases 
with reduced survey effort (Oyafuso et al., 2021) can per- 
haps be mitigated by the use of spatiotemporal model- 
based methods, such as those in the VAST package. Similar 
improvements were also observed when simple random 
sampling was used with the model-based estimator. 

In a related study, in which the effect of sampling den- 
sity changes on biomass estimates was examined for the 
same 3 species analyzed here, the model-based estimator 
when used with the stratified-random sampling design 
performed better than the design-based estimator for all 
3 species at 4 different levels of survey effort (von Szalay 
et al.”). In other studies, the advantages of a spatiotempo- 
ral model-based estimator have been demonstrated more 


2 von Szalay, P. G., E. A. Laman, S. Kotwicki, L. A. K. Barnett, and 
K. Ono. In preparation. Quantifying the effects of sample size 
and species distribution on the precision and accuracy of abun- 
dance estimates from bottom trawl surveys in the Gulf of Alaska. 


broadly. O’Leary et al. (2020) found that the spatial and 
temporal scope of total biomass for walleye pollock (G. chal- 
cogrammus) in the northern Bering Sea was improved in a 
stock assessment model when a model-based biomass index 
from multiple surveys was used as model input in combina- 
tion with the standard design-based estimates. In addition, 
stock assessment models for dusky rockfish (S. ciliatus) 
and northern rockfish (S. polyspinis) in the GOA have 
incorporated such model-based estimates because they are 
less noisy than the sampling-design-based CVs (Lunsford 
et al., 2015; Cunningham et al., 2018). 

The improvements we observed with the model-based 
estimator are likely due to an increase in the effective sam- 
ple size when data from multiple surveys are used to derive 
biomass estimates. Thorson and Haltuch (2019) showed 
that the effective sample size increased by 17% and that the 
model fit was better, resulting in smaller standard errors of 
estimated spawning biomass, when a spatiotemporal model- 
based estimator was used in a stock synthesis assessment 
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Figure 8 


Relative performance of 7 scenarios with different combinations of estimators and sampling 
designs in estimating the variance of biomass for Pacific ocean perch (Sebastes alutus) based on 
100 simulated bottom-trawl surveys in the Gulf of Alaska. Each box plot represents the distri- 
bution of a performance metric, (A) coefficient of variation, (B) relative bias, or (C) relative root- 
mean-square error (RMSE), for a scenario over the 10 survey years between 1996 and 2015. The 
model-based estimator, the vector autoregressive spatiotemporal (VAST) model in the R package 
VAST, was paired with 4 sampling designs: systematic (VSY scenario), simple random (VRS sce- 
nario), stratified random based on information scores (VIS scenario), and traditional stratified 
random (VTRS scenario). The designed-based estimator was paired with 3 sampling designs: sim- 
ple random (SRS scenario), stratified random based on information scores (WIS scenario), and tra- 
ditional stratified random based on Neyman allocation (TRS scenario). In each box plot, the thick 
black line is the median. The lower and upper hinges of the box plots correspond to the 25th and 
75th percentiles. The whiskers extend from the hinge to the farthest value no more than 1.5 times 
the interquartile range (distance between the 25th and 75th percentiles). Data points beyond the 


end of the whiskers, outliers, are plotted individually. 


model for lingcod (Ophiodon elongatus) in the California 
Current. However, we acknowledge that additional work 
with more than the 3 species examined in our study is 
needed before our results can be generalized. It is possible 
that the VAST estimator does not perform as well for spe- 
cies with a more limited spatial distribution. 

Furthermore, both bias and variance tended to be high 
for the model-based estimator when the sampling design 
is highly unbalanced. This was the case for all 3 species 
when the VAST estimation model was paired with the 
IS-based sampling design (VIS scenario), which concen- 
trated sampling density in areas of high biomass. 


Recommendation to develop stratification methods 
based on information scores 


For all 3 species, [S-based stratification in combination 
with design-based estimation yielded results that were by 
far better than any other method. The main reason this 


method outperformed the traditional design schemes (ran- 
dom and systematic) was that stratification was informed 
by historical data in the form of ISes, which can be described 
as estimated expectation of population abundance and its 
temporal variability. Of course, for this approach to be valid, 
one must be confident in the predictions of density based 
on historical data to adequately represent current and per- 
haps future species distributions. In contrast, the stratum 
boundaries in traditional random sampling are defined by 
depth and management area. Using ISes to stratify the 
survey area resulted in strata with more homogenous vari- 
ances than the strata in the traditional Neyman allocation. 
This difference is due to the variance component of the IS, 
which increases the efficiency of IS-based stratification 
compared to that of the traditional stratification scheme. 
Furthermore, in contrast to using Neyman allocation, the 
use of ISes ensures that samples are representative in 
areas where density is persistently high and temporal vari- 
ance is low. In Neyman allocation, strata with high density 
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but low variance are apportioned a relatively small number 
of stations. Such allocation is undesirable for fisheries sur- 
veys because, in addition to precise estimates of abundance, 
representative samples for other population measures, 
such as age and size compositions, are also needed. In the 
extreme case of the abundance of a stratum being very high 
but uniform, the Neyman strategy would allocate no sam- 
ples to that stratum. 

In the future, it may be problematic to use historical data 
to define stratum boundaries under the IS approach if fish 
distributions change in response to climate change as they 
are in many regions (Hobday and Evans, 2013; Maureaud 
et al., 2021; Hollowed et al., 2022). If the shift in distribu- 
tions is relatively gradual over time, it may still be appropri- 
ate to use relatively recent data while excluding older data 
to define stratum boundaries on the basis of ISes. However, 
if the shift in distribution patterns is more abrupt, it may 
not be appropriate to use any sampling design that defines 
strata on the basis of historical data, as is the case for both 
the WIS and TRS approaches. In this situation, it may be 
necessary to conduct a pilot survey prior to the main survey 
to establish new stratum boundaries. 

The finding that the design based on ISes outperformed 
the traditional sampling design is not surprising given 
that the stratification was performed individually for each 
species in the case of the WIS scenario. Therefore, this 
method would be applicable to only single-species surveys. 
In contrast, the stratification scheme used for the TRS 
scenario was first optimized for each species, but because 
the GOA survey is a multispecies survey, the realized sta- 
tion allocation to each stratum was based on a weighted 
average of several principal species. This weighting effec- 
tively resulted in less than optimal station allocation for 
each species. The comparison is nevertheless informative 
because it illustrates how much better the biomass and 
variance estimates are if we can focus on sampling indi- 
vidual species, given their importance to the fishery. Per- 
haps it has been a mistake to make so many compromises 
with the traditional weighting across too many species. 

More research is needed to investigate how the relative 
performance of the WIS scenario compares with the TRS 
scenario in the context of a multispecies survey, in which 
the IS-based station allocation scheme would be subjected 
to the same compromises among several principal species 
as the traditional stratified-random sampling design. We 
envision generating multispecies ISes, which are linear 
combinations of the single-species scores associated with 
each of the principal species included in the TRS scenrio. 
Each component of these combinations could be weighted 
by a number of different factors important to manage- 
ment (e.g., the commercial value of each species, ecologi- 
cal importance of each species, or non-biological elements, 
such as cultural or political considerations). 

Although both the absolute and relative performance 
of a multispecies version of the WIS would necessar- 
ily be lower than the performance of the single-species 
version, we expect it to perform better than the TRS 
scenario because the stratum boundaries and station allo- 
cation scheme are strictly driven by historical data and 
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determined by the objective factors of fish density and its 
variance. Furthermore, it is likely that the IS approach can 
be improved by optimizing for both the relative weighting 
of the variance component and the number of strata used. 
In contrast, in the TRS scenario, strata are defined by ad 
hoc boundaries and depth of the management area, likely 
reducing the efficiency of the traditional sampling design. 
The allocation scheme of the TRS scenario also includes a 
penalty against strata that are relatively costly to sample 
in terms of time. This cost penalty may be problematic for 
estimating biomass and its variance for slope species, such 
as Pacific ocean perch and sablefish (Anoplopoma fim- 
bria), which primarily inhabit these high-cost strata and 
are also more patchily distributed than species, such as 
Pacific cod and arrowtooth flounder, that inhabit the con- 
tinental shelf. As a compromise with the objective to pre- 
serve the TRS approach, it may be feasible in the future to 
introduce a hybrid sampling design in which sampling for 
some species is done with the traditional design and sam- 
pling for other species, such as Pacific ocean perch, may be 
enhanced by sampling more in areas with high [Ses. 

The relatively poor performance of the model-based VIS 
scenario can likely be attributed to the unbalanced sam- 
pling design resulting from stratification based on ISes. 
Under this stratification scheme, in which the distributions 
of the [Ses were highly skewed (Suppl. Fig. 1) (online only), 
the bulk of the stations were assigned to just a few strata 
while the minimum of 2 required stations were assigned to 
all of the remaining strata, concentrating the spatial distri- 
bution of samples used to extrapolate to the entire survey 
area. The performance of the VIS in estimating biomass was 
particularly poor for Pacific ocean perch because this spe- 
cies is primarily confined to a depth range of approximately 
150-300 m, which exists only in a relatively narrow band 
on the upper continental slope of the survey area. Pacific 
cod and arrowtooth flounder, on the other hand, are more 
evenly distributed over a larger portion of the survey area 
(primarily over the relatively wide continental shelf); there- 
fore, the effect of the unbalanced sampling design is less 
pronounced for them. The IS-based sampling design tailors 
the station allocation to individual strata, resulting in better 
performance than the traditional stratification scheme with 
designed-based estimators because stations are allocated 
to areas where fish are most likely to be, thereby reducing 
the variance in biomass estimates while having no effect on 
bias. However, it does not work well with the model-based 
estimator, which is sensitive to the highly unbalanced sam- 
pling dictated by the IS-based stratification. 


Conclusions 


The results of our analysis of 6 scenarios with different 
combinations of sampling designs and estimators, in com- 
parison to the method used for the bottom-trawl survey 
conducted in the GOA since 1984, not only identify which 
scenario consistently outperformed the traditional approach 
but also indicate that the traditional approach of using a 
stratified-random sampling design with a design-based 
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estimator is not optimal for any of the 3 species consid- 
ered in this study. Estimates of abundance can likely be 
improved immediately by using the model-based estimator 
in the VAST package in combination with the traditional 
stratified-random sampling design (VTRS scenario). In the 
long term, stratification based on the IS proposed in this 
study, when used with a design-based estimator, could be a 
useful alternative to the traditional approach in improving 
estimates of abundance for species, such as Pacific ocean 
perch, with spatially restrictive or patchy distributions. 
However, in subsequent work, interspecific tradeoffs in a 
multispecies context will be required to justify reconsid- 
ering the traditional sampling design. Nonetheless, this 
approach of combining an IS-based sampling design with a 
design-based Wakabayashi estimator is promising, particu- 
larly for single-species surveys or for multispecies surveys 
in areas with few interspecific differences in distribution. 


Resumen 


Desde 1984, se han utilizado prospecciones de arrastre de 
fondo con un diseno de muestreo aleatorio estratificado 
para la evaluacién de poblaciones de especies de impor- 
tancia comercial en el Golfo de Alaska. Se evalu6é un nuevo 
diseno de muestreo estratificado para determinar si su uso 
podria mejorar la precisién y exactitud de las estimaciones 
de abundancia. En el enfoque propuesto para definir los 
estratos, se utilizaron datos de estudios histéricos para gen- 
erar lo que denominamos puntajes de informacion (ISes). 
Para determinar si el enfoque existente es 6ptimo, se com- 
paro el esquema de estratificacién tradicional con el nuevo 
método y otros 2 disefios de muestreo, utilizando tanto un 
estimador basado en el diseno como un estimador basado en 
el modelo con cada diseno. La robustez estadistica, medida 
en términos del coeficiente de variacién, sesgo y la raiz del 
error cuadratico medio, se compar6 entre 7 escenarios con 
diferentes combinaciones de estimadores y disenos de mues- 
treo mediante simulacién con un modelo lineal general- 
izado mixto espaciotemporal condicionado a observaciones 
historicas de capturas por unidad de esfuerzo de 3 especies. 
La combinacion del estimador basado en el disefno con el 
esquema de estratificacién basado en los ISes fue el mejor 
escenario en todas las medidas de desempeno para todas 
las especies. Este escenario consistentemente presento la 
varianza mas baja y el menor error total, y en general, no 
estuvo sesgado. En contraste, la comparacion del estimador 
basado en el modelo con este disenio de muestreo fue por 
mucho, el escenario con peor desempeno. El desempeno del 
enfoque existente fue promedio. 
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and Ronald Goldberg 
Temperature-related changes in species composition of 
juvenile finfish on a rock reef in Long Island Sound 
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Corrections 


On Page 280, Figure 4 should read as follows, with corrected 
values on the y-axis of graphs in all panels except in the top 
2 panels for black sea bass and cunner: 
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Figure 4 


Abundance of juveniles for dominant finfish species caught in traps on a cobble and boulder reef in Long 
Island Sound during June, July, and August in 2004—2008 and in 2016. Catch per unit of effort (CPUE), 
used as an index of abundance, was calculated by summing the number of all fish caught and dividing by 
the total annual fishing effort (number of traps multiplied by number of sampling trips). The dominant 
species were the black sea bass (Centropristis striata), cunner (Tautoga adspersus), tautog (Tautoga oni- 
tis), rock gunnel (Pholis gunnellus), grubby (Myoxocephalus aenaeus), scup (Stenotomus chrysops), and 
oyster toadfish (Opsanus tau). Note that the graphs for black sea bass and cunner have a different range 
of CPUE values on the y-axis (up to 0.6) than graphs for the other species (values up to 0.06). 
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followed by initials for first name and, if given, mid- 
dle name of first author; then list names of additional 
authors with initials before last names). Year. Title of 
article. Abbreviated title of the journal in which it was 
published. Always include either the range of page num- 
bers (for a journal article) or a total number of pages (for 
a book or other type of publication). List a sequence of 
citations in the general text chronologically, for example, 
“(Smith, 1932; Green, 1998; Smith and Jones, 2015).” 


Acknowledgments should be no more than 6 lines of 
text. Only those who have contributed in an outstanding 
way should be acknowledged by name. For recognition of 
other persons or groups, use a general term, such as crew, 
observers, or research coordinators, and do not include 
names with these terms. 


Digital object identifier (doi) code ensures that a publica- 
tion has a permanent location online. A doi link (which 
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may include a doi code) should be included at the end of 
citations of published literature. Authors are responsi- 
ble for submitting accurate doi links. Faulty links will be 
deleted at the page-proof stage. 


Footnotes are used for all documents that have not been 
formally peer reviewed and for observations and personal 
communications, but these types of references should be 
cited sparingly in manuscripts submitted to the journal. 

All reference documents, administrative reports, inter- 
nal reports, progress reports, project reports, contract 
reports, personal observations, personal communications, 
unpublished data, manuscripts in review, and council meet- 
ing notes are footnoted in 10-point font and placed at the 
bottom of the page on which they are first cited. Footnote 
format is the same as that for formal literature citations. A 
link to the online source (e.g., [Available from http://www... , 
accessed July 2017.]), or the mailing address of the agency 
or department holding the document, should be provided 
so that readers may obtain a copy of the document. 


Tables are often overused in scientific papers; it is seldom 
necessary to present all the data associated with a study. 
Tables should not be excessive in size and must be cited 
in numerical order in the text. Headings should be short 
but ample enough to allow the table to be intelligible on 
its own. 

All abbreviations and unusual symbols must be 
explained in the table legend. Other incidental com- 
ments may be footnoted with numeral footnote markers. 
Use asterisks only to indicate significance in statistical 
data. Do not put a table legend on a page separate from 
the table; place the legend above the table. Do not submit 
tables in photo mode. 


e Note probability with a capital, italic P. 


e Provide a zero before all decimal points for values less 
than one (e.g., 0.07). 


e Round all values to 2 decimal points. 


e Use a comma in numbers of 5 digits or more (e.g., 
13,000 but 3000). 


Figures must be cited in numerical order in the text. 
Graphics should aid in the comprehension of the text, but 
they should be limited to presenting patterns rather than 
raw data. The number of figures should not exceed 1 figure 
for every 4 pages of text. 

Figure legends should explain all symbols and abbrevi- 
ations seen in the figure and should be double spaced on a 
separate page at the end of the manuscript. 

Line art and halftone figures should be saved at res- 
olutions >600 dots per inch (dpi) and >300 dpi, respec- 
tively. Color is allowed in figures to show morphological 
differences among species (for species identification), to 
show stain reactions, to show gradations (such as those of 
temperature and salinity within maps), and to distinguish 
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between numerous lines and symbols in graphs. Figures 
approved for color should be saved in CMYK format. 

All figures must be submitted as PDF, TIFF, or EPS 
files. 


e Capitalize only the first letter of the first word and 
proper nouns in all labels within figures. 


¢ Do not use overly large font sizes for labels in maps and 
for axis labels in graphs. 


e Use the same point size for all labels, except for panel 
labels (e.g., A), which should be slightly larger than 
other labels (e.g., 11 point versus 8 point). 


e Use asans serif font for all labels. Panel labels (e.g., A), 
however, should be in Times New Roman font. 


e Do not use bold fonts or bold lines in figures. Do not 
use italic fonts. Exceptions include use of italic fonts 
for labels of bodies of water in maps and a bold font for 
panel labels (e.g., A). 


e Do not place outline rules around graphs. 


e Do not include vertical and horizontal lines in the back- 
ground of graphs. Ticks for values on the x-axis and 
y-axis will suffice. 


e Place a north arrow and label degrees latitude and lon- 
gitude (e.g., 170°E) in all maps. If scale of map requires 
more than degrees, use degrees minutes, not decimal 
degrees. 


e Place panel labels (e.g., A, B, C) within the upper-left 
area of each graph or photo in a multi-panel figure, 
from left to right, then top to bottom. If the letter is 
not visible against a dark background, put a white box 
behind it. Do not use white labels. 


e Avoid placing labels vertically or diagonally. Y-axis 
labels can be vertical. Words in horizontal labels can be 
stacked vertically to fit. 


e Use symbols, shading, or patterns (not clip art) in maps 
and graphs. 


e For scale bars in maps, use kilometers. Use the label 
km or kilometers (lowercased). 


Supplementary materials that are considered essential, 
but are too large or impractical for inclusion in a paper 
(e.g., metadata, figures, tables, videos, and websites), may 
be provided at the end of an article. These materials are 
subject to the editorial standards of the journal. A URL to 
the supplementary material and a brief explanation for 
including such material should be sent at the time of ini- 
tial submission of the paper to the journal. 


e Metadata, figures, and tables should be submitted in 
standard digital format (MS Word or PDF file) and 
should be cited in the general text, for example, as 
“.. was determined (Suppl. Table 3, Suppl. Fig. 1).” 


e Websites should be cited with a URL in the general 
Text: 


e Videos must not be larger than 30 MB to allow a swift 
technical response for viewing the video. Authors 
should consider whether a short video uniquely cap- 
tures what text alone cannot capture for the under- 
standing of a process or behavior under examination 
in the article. Supply an online link to the location of 
the video. 


Copyright law does not apply to Fishery Bulletin, which 
falls within the public domain. However, if an author 
reproduces any part of an article from Fishery Bulletin, 
reference to source is considered correct form (e.g., Source: 
Fish. Bull. 117:105). 


Failure to follow these guidelines 
and failure to correspond with editors 
in a timely manner will delay 
publication of a manuscript. 


Submission of manuscript 


Submit a manuscript online at the ScholarOne Manu- 
scripts website for Fishery Bulletin. Commerce Depart- 
ment authors must provide proof of internal clearance 
of their manuscripts with either a completed and signed 
NOAA Form 25-700 or a copy of the clearance email from 
the Research Publication Tracking System. For further 
details on electronic submission, please contact the asso- 
ciate editor, Cara Mayo, at 


cara.mayo@noaa.gov. 


When requested, the text and tables should be submitted 
in MS Word format. Each figure should be sent as a sep- 
arate PDF, TIFF, or EPS file. Send a copy of a figure in 
the original software if conversion to any of these formats 
yields a degraded version of the figure. 


Questions? If you have questions regarding these guide- 
lines, please contact the managing editor, Kathryn 
Dennis, at 


kathryn.dennis@noaa.gov. 


Questions regarding manuscripts under review should be 
addressed to Associate Editor Cara Mayo. 
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